require(estimatr)
require(doParallel)

#### data preprocessing ####
## Korean data
Korea.data <- read.csv("Korea_data.csv")
Korea.data$C  <- (Korea.data$Gubun == 3) * 1
Korea.data$T1 <- (Korea.data$Gubun == 1) * 1
Korea.data$T2 <- (Korea.data$Gubun == 2) * 1
Korea.data$condition <- ifelse(Korea.data$Gubun == 3, 1,
                               ifelse(Korea.data$Gubun == 1, 2, 3))
Korea.data$male <- ifelse(Korea.data$Q4 == 1, 1,
                          ifelse(Korea.data$Q4 == 2, 0, NA))
Korea.data$gender <- Korea.data$male + 1
Korea.data$age <- Korea.data$Q2
Korea.data$educ <- Korea.data$Q5
Korea.data$educ.cat <- ifelse(Korea.data$Q5 < 3, 1,
                              ifelse(Korea.data$Q5 == 3, 2,
                                     ifelse(Korea.data$Q5 > 3, 3, NA)))
Korea.data$middle.educ <- (Korea.data$educ.cat == 2) * 1
Korea.data$high.educ <- (Korea.data$educ.cat == 3) * 1
Korea.data$income <- ifelse(Korea.data$Q9 == 11, NA, Korea.data$Q9)
Korea.data$high.income <- ifelse(Korea.data$Q9 == 11, NA,
                                 (Korea.data$Q9 >= median(Korea.data$Q9[Korea.data$Q9 != 11])) * 1)
Korea.data$income.cat <- cut(Korea.data$income, 4, labels = 1:4)
Korea.data$income.incl.NA <- factor(ifelse(is.na(Korea.data$income.cat), 5,
                                           Korea.data$income.cat))
Korea.data$ideology <- ifelse(Korea.data$Q11 < 3, 5,
                              ifelse(Korea.data$Q11 == 3 | Korea.data$Q11 == 4, 4,
                                     ifelse(Korea.data$Q11 == 5, 3,
                                            ifelse(Korea.data$Q11 == 6 | Korea.data$Q11 == 7, 2,
                                                   ifelse(Korea.data$Q11 > 7, 1, NA)))))
Korea.data$BS <- rowMeans(cbind(Korea.data$Q7_1, Korea.data$Q7_2, 
                                Korea.data$Q7_3))
Korea.data$HS <- rowMeans(cbind(Korea.data$Q7_4, Korea.data$Q7_5, 
                                Korea.data$Q7_6))
Korea.data$MS <- rowMeans(cbind(Korea.data$Q7_7, 6 - Korea.data$Q7_8, 
                                Korea.data$Q7_9))
Korea.data$gender.gap.est <- Korea.data$Q16_1
Korea.data$gender.rank.est <- Korea.data$Q16_2
Korea.data$income.gap.est <- Korea.data$Q17_1
Korea.data$income.rank.est <- Korea.data$Q17_2
Korea.data$gender.gap.undest <- (Korea.data$Q16_1 > 68.8) * 1
Korea.data$income.gap.undest <- (Korea.data$Q17_1 > 20) * 1
Korea.data$gender.rank.undest <- (Korea.data$Q16_2 > 1) * 1
Korea.data$income.rank.undest <- (Korea.data$Q17_2 > 1) * 1
Korea.data$gender.gap.overest <- (Korea.data$Q16_1 < 68.8) * 1
Korea.data$income.gap.overest <- (Korea.data$Q17_1 < 20) * 1
Korea.data$gender.rank.overest <- (Korea.data$Q16_2 == 1) * 1
Korea.data$income.rank.overest <- (Korea.data$Q17_2 == 1) * 1
Korea.data$intv.gender.gap <- Korea.data$Q21_3
Korea.data$intv.gender.policy <- Korea.data$Q21_1
Korea.data$intv.income.gap <- Korea.data$Q21_4
Korea.data$intv.income.policy <- Korea.data$Q21_2
Korea.data$feeling.women <- Korea.data$Q22_1
Korea.data$feeling.men <- Korea.data$Q22_2
Korea.data$feeling.gender.diff <- Korea.data$Q22_1 - Korea.data$Q22_2
Korea.data$feeling.poor <- Korea.data$Q22_3
Korea.data$feeling.rich <- Korea.data$Q22_4
Korea.data$feeling.income.diff <- Korea.data$Q22_3 - Korea.data$Q22_4

## Japanese data
Japan.data <- read.csv("Japan_data.csv")
Japan.data$C <- (Japan.data$condition == 1) * 1
Japan.data$T1 <- (Japan.data$condition == 2) * 1
Japan.data$T2 <- (Japan.data$condition == 3) * 1
Japan.data$age <- Japan.data$Q3.1 + 16
Japan.data$male <- ifelse(Japan.data$Q4.3 == 3, NA,
                          (Japan.data$Q4.3 == 2) * 1)
Japan.data$gender <- Japan.data$Q4.3
Japan.data$educ <- Japan.data$Q4.5
Japan.data$educ.cat <- ifelse(Japan.data$Q4.5 < 3, 1, 
                              ifelse(Japan.data$Q4.5 < 6, 2, 3))
Japan.data$middle.educ <- (Japan.data$educ.cat == 2) * 1
Japan.data$high.educ <- (Japan.data$educ.cat == 3) * 1
Japan.data$income <- ifelse(Japan.data$Q6.3 == 11, NA, Japan.data$Q6.3)
Japan.data$high.income <- ifelse(Japan.data$Q6.3 == 11, NA, 
                                 (Japan.data$Q6.3 >= median(Japan.data$income, na.rm = TRUE)) * 1)
Japan.data$income.cat <- cut(Japan.data$income, 4, labels = 1:4)
Japan.data$income.incl.NA <- factor(ifelse(is.na(Japan.data$income.cat), 5, 
                                           Japan.data$income.cat))
Japan.data$ideology <- ifelse(is.na(Japan.data$Q8.1), 3, 
                              ifelse(Japan.data$Q8.1 < 3, 5, 
                                     ifelse(Japan.data$Q8.1 < 5, 4, 
                                            ifelse(Japan.data$Q8.1 == 5, 3, 
                                                   ifelse(Japan.data$Q8.1 < 8, 2, 1)))))
Japan.data$BS <- rowMeans(cbind(Japan.data$Q5.1_1, Japan.data$Q5.1_2, 
                                Japan.data$Q5.1_3))
Japan.data$HS <- rowMeans(cbind(Japan.data$Q5.1_4, Japan.data$Q5.1_5, 
                                Japan.data$Q5.1_6))
Japan.data$MS <- rowMeans(cbind(Japan.data$Q5.1_7, 6 - Japan.data$Q5.1_8, 
                                Japan.data$Q5.1_9))
Japan.data$gender.gap.est <- Japan.data$Q11.1_1
Japan.data$gender.rank.est <- 5 - Japan.data$Q11.3
Japan.data$income.gap.est <- Japan.data$Q12.1_1
Japan.data$income.rank.est <- 5 - Japan.data$Q12.3
Japan.data$gender.gap.undest <- (Japan.data$gender.gap.est >= 79) * 1
Japan.data$gender.rank.undest <- (Japan.data$gender.rank.est > 1) * 1
Japan.data$income.gap.undest <- (Japan.data$income.gap.est >= 20) * 1
Japan.data$income.rank.undest <- (Japan.data$income.rank.est > 1) * 1
Japan.data$gender.gap.overest <- 1 - Japan.data$gender.gap.undest
Japan.data$gender.rank.overest <- 1 - Japan.data$gender.rank.undest
Japan.data$income.gap.overest <- 1 - Japan.data$income.gap.undest
Japan.data$income.rank.overest <- 1 - Japan.data$income.rank.undest
Japan.data$intv.gender.gap <- Japan.data$Q15.1_3
Japan.data$intv.gender.policy <- Japan.data$Q15.1_1
Japan.data$intv.income.gap <- Japan.data$Q15.1_4
Japan.data$intv.income.policy <- Japan.data$Q15.1_2
Japan.data$feeling.women <- Japan.data$Q15.3_1
Japan.data$feeling.men <- Japan.data$Q15.3_2
Japan.data$feeling.gender.diff <- Japan.data$feeling.women - Japan.data$feeling.men
Japan.data$feeling.poor <- Japan.data$Q15.3_3
Japan.data$feeling.rich <- Japan.data$Q15.3_4
Japan.data$feeling.income.diff <- Japan.data$feeling.poor - Japan.data$feeling.rich
Japan.data$block <- Japan.data$Q4.3 * 100 + 
  Japan.data$gender.gap.undest * 10 + Japan.data$income.gap.undest

combined.data <- rbind(cbind(Japan = 0, 
                             Korea.data[, c("T1", "T2", "age", "gender", "male", 
                                            "educ.cat", "income.incl.NA", 
                                            "ideology", "HS", "BS", "MS", 
                                            "gender.gap.undest", "income.gap.undest", 
                                            "gender.rank.undest", "income.rank.undest", 
                                            "intv.gender.gap", "intv.gender.policy", 
                                            "feeling.women", "feeling.gender.diff", 
                                            "intv.income.gap", "intv.income.policy",  
                                            "feeling.poor", "feeling.income.diff")]), 
                       cbind(Japan = 1, 
                             Japan.data[, c("T1", "T2", "age", "gender", "male", 
                                            "educ.cat", "income.incl.NA", 
                                            "ideology", "HS", "BS", "MS", 
                                            "gender.gap.undest", "income.gap.undest", 
                                            "gender.rank.undest", "income.rank.undest", 
                                            "intv.gender.gap", "intv.gender.policy", 
                                            "feeling.women", "feeling.gender.diff", 
                                            "intv.income.gap", "intv.income.policy",  
                                            "feeling.poor", "feeling.income.diff")]))

#### user-defined functions ####
# effect.plot: plot point estimates and 95% CIs from lm_robust models
# Note: When group is specified, compute subgroup CATEs and 
# differences using linear combinations of coefficients.
effect.plot <- function(data, outcome, treatment, group = NULL, 
                        FE, ylim, yaxis.labs, title, 
                        group.labs = NULL, 
                        cex.1 = 0.8, cex.2 = 0.6) {
  plot(NULL, NULL, bty = "n", xlim = c(0.5, 3.5), ylim = ylim, 
       xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  abline(h = 0, col = "gray80")
  abline(h = yaxis.labs[yaxis.labs != 0], lty = 3, col = "gray80")
  axis(2, at = yaxis.labs, lwd = 0.5)
  if (is.null(group)) {
    model <- formula(paste(outcome, "~ T1 + T2"))
    if (FE) {
      lm.out <- lm_robust(model, data, fixed_effects = block)
    } else {
      lm.out <- lm_robust(model, data)
    }
    segments(1.5, lm.out$conf.low["T1"], 1.5, lm.out$conf.high["T1"])
    points(1.5, lm.out$coefficients["T1"], 
           pch = ifelse(treatment == "T1", 19, 21), 
           bg = ifelse(treatment == "T1", NA, "white"))
    segments(2.5, lm.out$conf.low["T2"], 2.5, lm.out$conf.high["T2"])
    points(2.5, lm.out$coefficients["T2"], 
           pch = ifelse(treatment == "T1", 21, 19), 
           bg = ifelse(treatment == "T1", "white", NA))
    mtext(title, line = 1, font = 2, cex = 0.8)
    mtext(c("T1", "T2"), side = 1, line = 0, at = c(1.5, 2.5), cex = 0.6)
  } else {
    model <- formula(paste(outcome, "~", treatment, "*", group))
    if (FE) {
      lm.out <- lm_robust(model, data, fixed_effects = block)
    } else {
      lm.out <- lm_robust(model, data)
    }
    CATE.est <- lm.out$coefficients[treatment] + 
      lm.out$coefficients[paste0(treatment, ":", group)]
    CATE.se <- sqrt(lm.out$vcov[treatment, treatment] + 
                      lm.out$vcov[paste0(treatment, ":", group), paste0(treatment, ":", group)] + 
                      2 * lm.out$vcov[treatment, paste0(treatment, ":", group)])
    CATE.ci.low <- CATE.est + qt(0.025, lm.out$df[1]) * CATE.se
    CATE.ci.high <- CATE.est + qt(0.975, lm.out$df[1]) * CATE.se
    segments(1, lm.out$conf.low[treatment], 1, lm.out$conf.high[treatment])
    points(1, lm.out$coefficients[treatment], pch = 19)
    segments(2, CATE.ci.low, 2, CATE.ci.high)
    points(2, CATE.est, pch = 19)
    segments(3, lm.out$conf.low[paste0(treatment, ":", group)], 
             3, lm.out$conf.high[paste0(treatment, ":", group)])
    points(3, lm.out$coefficients[paste0(treatment, ":", group)], pch = 4)
    mtext(title, line = 0, font = 2, cex = cex.1)
    mtext(c(group.labs, "diff"), side = 1, line = 0, at = 1:3, cex = cex.2)
  }
}

## effect.plot.2: plot subgroup-specific treatment effects using interaction terms
effect.plot.2 <- function(data, outcome, treatment, group.1, group.2, 
                          ylim, yaxis.labs, title, 
                          group.labs) {
  plot(NULL, NULL, bty = "n", xlim = c(0.5, 4.5), ylim = ylim, 
       xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  abline(h = 0, col = "gray80")
  abline(h = yaxis.labs[yaxis.labs != 0], lty = 3, col = "gray80")
  axis(2, at = yaxis.labs, lwd = 0.5)
  model <- formula(paste(outcome, "~", treatment, "*", group.1, "*", group.2))
  lm.out <- lm_robust(model, data)
  var.1 <- paste0(treatment, ":", group.1)
  var.2 <- paste0(treatment, ":", group.2)
  var.3 <- paste0(treatment, ":", group.1, ":", group.2)
  CATE.est.0.0 <- lm.out$coefficients[treatment]
  CATE.ci.low.0.0 <- lm.out$conf.low[treatment]
  CATE.ci.high.0.0 <- lm.out$conf.high[treatment]
  CATE.est.1.0 <- lm.out$coefficients[treatment] + 
    lm.out$coefficients[var.1]
  CATE.se.1.0 <- sqrt(lm.out$vcov[treatment, treatment] + 
                        lm.out$vcov[var.1, var.1] + 
                        2 * lm.out$vcov[treatment, var.1])
  CATE.ci.low.1.0 <- CATE.est.1.0 + qt(0.025, lm.out$df[1]) * CATE.se.1.0
  CATE.ci.high.1.0 <- CATE.est.1.0 + qt(0.975, lm.out$df[1]) * CATE.se.1.0
  CATE.est.0.1 <- lm.out$coefficients[treatment] + 
    lm.out$coefficients[var.2]
  CATE.se.0.1 <- sqrt(lm.out$vcov[treatment, treatment] + 
                        lm.out$vcov[var.2, var.2] + 
                        2 * lm.out$vcov[treatment, var.2])
  CATE.ci.low.0.1 <- CATE.est.0.1 + qt(0.025, lm.out$df[1]) * CATE.se.0.1
  CATE.ci.high.0.1 <- CATE.est.0.1 + qt(0.975, lm.out$df[1]) * CATE.se.0.1
  CATE.est.1.1 <- lm.out$coefficients[treatment] + 
    lm.out$coefficients[var.1] + 
    lm.out$coefficients[var.2] + 
    lm.out$coefficients[var.3]
  CATE.se.1.1 <- sqrt(lm.out$vcov[treatment, treatment] + 
                        lm.out$vcov[var.1, var.1] + 
                        lm.out$vcov[var.2, var.2] + 
                        lm.out$vcov[var.3, var.3] + 
                        2 * lm.out$vcov[treatment, var.1] + 
                        2 * lm.out$vcov[treatment, var.2] + 
                        2 * lm.out$vcov[treatment, var.3] + 
                        2 * lm.out$vcov[var.1, var.2] + 
                        2 * lm.out$vcov[var.1, var.3] + 
                        2 * lm.out$vcov[var.2, var.3])
  CATE.ci.low.1.1 <- CATE.est.1.1 + qt(0.025, lm.out$df[1]) * CATE.se.1.1
  CATE.ci.high.1.1 <- CATE.est.1.1 + qt(0.975, lm.out$df[1]) * CATE.se.1.1
  diff.est.0 <- lm.out$coefficients[var.1]
  diff.ci.low.0 <- lm.out$conf.low[var.1]
  diff.ci.high.0 <- lm.out$conf.high[var.1]
  diff.est.1 <- lm.out$coefficients[var.1] + 
    lm.out$coefficients[var.3]
  diff.se.1 <- sqrt(lm.out$vcov[var.1, var.1] + 
                      lm.out$vcov[var.3, var.3] + 
                      2 * lm.out$vcov[var.1, var.3])
  diff.ci.low.1 <- diff.est.1 + qt(0.025, lm.out$df[1]) * diff.se.1
  diff.ci.high.1 <- diff.est.1 + qt(0.975, lm.out$df[1]) * diff.se.1
  diff.in.diffs.est <- lm.out$coefficients[var.3]
  diff.in.diffs.low <- lm.out$conf.low[var.3]
  diff.in.diffs.high <- lm.out$conf.high[var.3]
  segments(0.85, CATE.ci.low.0.0, 0.85, CATE.ci.high.0.0)
  points(0.85, CATE.est.0.0, pch = 19)
  segments(1.15, CATE.ci.low.0.1, 1.15, CATE.ci.high.0.1)
  points(1.15, CATE.est.0.1, pch = 15)
  segments(1.85, CATE.ci.low.1.0, 1.85, CATE.ci.high.1.0)
  points(1.85, CATE.est.1.0, pch = 19)
  segments(2.15, CATE.ci.low.1.1, 2.15, CATE.ci.high.1.1)
  points(2.15, CATE.est.1.1, pch = 15)
  segments(2.85, diff.ci.low.0, 2.85, diff.ci.high.0)
  points(2.85, diff.est.0, pch = 21, bg = "white")
  segments(3.15, diff.ci.low.1, 3.15, diff.ci.high.1)
  points(3.15, diff.est.1, pch = 22, bg = "white")
  segments(4, diff.in.diffs.low, 4, diff.in.diffs.high)
  points(4, diff.in.diffs.est, pch = 4)
  mtext(title, line = 0, font = 2, cex = 0.8)
  mtext(c(group.labs, "diff b/w\ngenders", "diff in\ndiffs"), 
        side = 1, line = 1, at = 1:4, cex = 0.6)
}

# equivalence.F.test: equivalence test for mean balance across multiple groups 
# using a noncentral F distribution
equivalence.F.test <- function (data,  # data frame containing variables
                                variable,  # variable to be tested
                                group,  # group variable
                                epsilon  # equivalence limits
) {
  variable.comp <- na.omit(data[, variable])
  group.comp <- data[! is.na(data[, variable]), group]
  N <- length(variable.comp)
  k <- length(unique(group.comp))
  n <- table(group.comp)
  n.bar <- mean(n)
  grand.mean <- mean(variable.comp)
  group.mean <- tapply(variable.comp, group.comp, mean)
  SS.between <- 0
  SS.within <- rep(NA, k)
  for (i in 1:k) {
    SS.between <- SS.between + 
      (n[i] / n.bar) * (group.mean[i] - grand.mean) ^ 2
    SS.within[i] <- sum((variable.comp[group.comp == i] - group.mean[i]) ^ 2)
  }
  psi.sq.hat <- SS.between / 
    (sum(SS.within) / (N - k))
  noncentrality.par <- n.bar * epsilon ^ 2
  p.value <- pf(psi.sq.hat / ((k - 1) / n.bar), k - 1, N - k, noncentrality.par)
  result <- data.frame(matrix(round(group.mean, 3), 1, k), p = round(p.value, 3))
  rownames(result) <- variable
  result
}

# t.test.fun: wrapper for two-sample t-tests that returns group means 
# and a p-value
t.test.fun <- function(formula, data) {
  results <- t.test(formula, data)
  out <- c(results$estimate, results$p.value)
  names(out) <- c("Group 0", "Group 1", "p value")
  round(out, 3)
}

# coef.table: format lm_robust output and compute "beta + delta" 
# (main effect + interaction term) for readability
coef.table <- function(x, blocking, T1) {
  out <- summary(x)$coefficients
  vcv <- x$vcov
  df <- x$df.residual
  n.vars <- nrow(out)
  vcv.rank <- nrow(vcv)
  int <- 1 - blocking * 1
  if (T1) {
    beta.plus.delta.est <- out[1 + int, 1] + out[n.vars - 1, 1]
    beta.plus.delta.se <- sqrt(vcv[1 + int, 1 + int] + vcv[vcv.rank - 1, vcv.rank - 1] + 2 * vcv[1 + int, vcv.rank - 1])
    beta.plus.delta.p <- 2 * (1 - pt(abs(beta.plus.delta.est) / beta.plus.delta.se, df))
  } else {
    beta.plus.delta.est <- out[2 + int, 1] + out[n.vars, 1]
    beta.plus.delta.se <- sqrt(vcv[2 + int, 2 + int] + vcv[vcv.rank, vcv.rank] + 2 * vcv[2 + int, vcv.rank])
    beta.plus.delta.p <- 2 * (1 - pt(abs(beta.plus.delta.est) / beta.plus.delta.se, df))
  }
  out <- rbind(out[c((1 + int):(3 + int), n.vars - 1, n.vars), c(1, 2, 4)], 
               c(beta.plus.delta.est, beta.plus.delta.se, beta.plus.delta.p), 
               c(x$nobs, NA, NA))
  rownames(out)[nrow(out) - 1] <- "beta + delta"
  rownames(out)[nrow(out)] <- "N"
  round(out, 3)
}

# two.way.anova.fun: helper for two-way ANOVA summaries (cell means, 
# within-group differences, and the interaction p-value)
two.way.anova.fun <- function(y, x1, x2, data) {
  results <- aov(formula(paste(y, "~", x1, "*", x2)), data)
  out <- c(mean(unlist(data[data[, x1] == 0 & data[, x2] == 0, y]), na.rm = TRUE), 
           mean(unlist(data[data[, x1] == 0 & data[, x2] == 1, y]), na.rm = TRUE), 
           mean(unlist(data[data[, x1] == 0 & data[, x2] == 1, y]), na.rm = TRUE) - 
             mean(unlist(data[data[, x1] == 0 & data[, x2] == 0, y]), na.rm = TRUE), 
           mean(unlist(data[data[, x1] == 1 & data[, x2] == 0, y]), na.rm = TRUE), 
           mean(unlist(data[data[, x1] == 1 & data[, x2] == 1, y]), na.rm = TRUE), 
           mean(unlist(data[data[, x1] == 1 & data[, x2] == 1, y]), na.rm = TRUE) - 
             mean(unlist(data[data[, x1] == 1 & data[, x2] == 0, y]), na.rm = TRUE), 
           summary(results)[[1]]$'Pr(>F)'[3])
  names(out) <- c("Group 0-0", "Group 0-1", "Diff in Group 0", 
                  "Group 1-0", "Group 1-1", "Diff in Group 1", "p value")
  round(out, 3)
}

#### distribution of the outcome variables ####
# Figure 1
cairo_pdf("Figure_1a.pdf", width = 3, height = 1.5, pointsize = 5)
par(mar = c(2, 4, 1, 1), lwd = 0.5)
hist(Korea.data$income.gap.est, main = "", xlab = "", axes = FALSE)
abline(v = 20)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
mtext("Actual value", line = 0, at = 20)
dev.off()

cairo_pdf("Figure_1b.pdf", width = 3, height = 1.5, pointsize = 5)
par(mar = c(2, 4, 1, 1), lwd = 0.5)
hist(Korea.data$income.rank.est, breaks = seq(0.5, 4.5, 1), 
     main = "", xlab = "", axes = FALSE)
axis(1, 1:4, c("1st\u201310th", "11th\u201320th", "21st\u201330th", "31st\u201338th"), 
     lwd = 0.5)
axis(2, lwd = 0.5)
mtext("Actual rank", line = 0, at = 1)
dev.off()

cairo_pdf("Figure_1c.pdf", width = 3, height = 1.5, pointsize = 5)
par(mar = c(2, 4, 1, 1), lwd = 0.5)
hist(Korea.data$gender.gap.est, main = "", xlab = "", axes = FALSE)
abline(v = 68.8)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
mtext("Actual value", line = 0, at = 68.8)
dev.off()

cairo_pdf("Figure_1d.pdf", width = 3, height = 1.5, pointsize = 5)
par(mar = c(2, 4, 1, 1), lwd = 0.5)
hist(Korea.data$gender.rank.est, breaks = seq(0.5, 4.5, 1), 
     main = "", xlab = "", axes = FALSE)
axis(1, 1:4, c("1st\u201310th", "11th\u201320th", "21st\u201330th", "31st\u201338th"), 
     lwd = 0.5)
axis(2, lwd = 0.5)
mtext("Actual rank", line = 0, at = 1)
dev.off()

# Figure 2
cairo_pdf("Figure_2a.pdf", width = 3, height = 1.5, pointsize = 5)
par(mar = c(2, 4, 1, 1), lwd = 0.5)
hist(Japan.data$income.gap.est, main = "", xlab = "", axes = FALSE)
abline(v = 19.2)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
mtext("Actual value", line = 0, at = 19.2)
dev.off()

cairo_pdf("Figure_2b.pdf", width = 3, height = 1.5, pointsize = 5)
par(mar = c(2, 4, 1, 1), lwd = 0.5)
hist(Japan.data$income.rank.est, breaks = seq(0.5, 4.5, 1), 
     main = "", xlab = "", axes = FALSE)
axis(1, 1:4, c("1st\u201310th", "11th\u201320th", "21st\u201330th", "31st\u201338th"), 
     lwd = 0.5)
axis(2, lwd = 0.5)
mtext("Actual rank", line = 0, at = 1)
dev.off()

cairo_pdf("Figure_2c.pdf", width = 3, height = 1.5, pointsize = 5)
par(mar = c(2, 4, 1, 0.5), lwd = 0.5)
hist(Japan.data$gender.gap.est, main = "", xlab = "", axes = FALSE)
abline(v = 78.7)
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
mtext("Actual value", line = 0, at = 78.7)
dev.off()

cairo_pdf("Figure_2d.pdf", width = 3, height = 1.5, pointsize = 5)
par(mar = c(2, 4, 1, 1), lwd = 0.5)
hist(Japan.data$gender.rank.est, breaks = seq(0.5, 4.5, 1), 
     main = "", xlab = "", axes = FALSE)
axis(1, 1:4, c("1st\u201310th", "11th\u201320th", "21st\u201330th", "31st\u201338th"), 
     lwd = 0.5)
axis(2, lwd = 0.5)
mtext("Actual rank", line = 0, at = 1)
dev.off()

#### treatment effects on the primary outcomes ####
# Figure 3
cairo_pdf("Figure_3_1.pdf", width = 6, height = 1.5, pointsize = 7)
par(mar = c(1, 2.5, 2, 0.5), lwd = 0.5, oma = c(0, 0, 2, 0))
layout(matrix(1:5, 1, 5, byrow = TRUE))
effect.plot(Korea.data, "intv.income.gap", "T2", 
            FE = FALSE, ylim = c(-0.5, 1), yaxis.labs = seq(-0.5, 1, 0.25), 
            title = "Entire sample")
effect.plot(subset(Korea.data, T1 == 0), 
            "intv.income.gap", "T2", "income.gap.undest", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by estimate\n(ratio)", 
            c("Overest.", "Underest."))
effect.plot(subset(Korea.data, T1 == 0), 
            "intv.income.gap", "T2", "income.rank.undest", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by estimate\n(rank)", 
            c("Accurate", "Underest."))
effect.plot(subset(Korea.data, T1 == 0 & income.gap.undest == 1), 
            "intv.income.gap", "T2", "high.income", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by income\n(Only underest. for ratio)", 
            c("Low", "High"))
effect.plot(subset(Korea.data, T1 == 0 & income.rank.undest == 1), 
            "intv.income.gap", "T2", "high.income", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by income\n(Only underest. for rank)", 
            c("Low", "High"))
title("Income gap intervention (Korea)", 
      cex.main = 1.5, line = 1, outer = TRUE)
dev.off()

cairo_pdf("Figure_3_2.pdf", width = 6, height = 1.5, pointsize = 7)
par(mar = c(1, 2.5, 2, 0.5), lwd = 0.5, oma = c(0, 0, 2, 0))
layout(matrix(1:5, 1, 5, byrow = TRUE))
effect.plot(Japan.data, "intv.income.gap", "T2", 
            FE = TRUE, ylim = c(-0.5, 1), yaxis.labs = seq(-0.5, 1, 0.25), 
            title = "Entire sample")
effect.plot(subset(Japan.data, T1 == 0), 
            "intv.income.gap", "T2", "income.gap.undest", 
            TRUE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by estimate\n(ratio)", 
            c("Overest.", "Underest."))
effect.plot(subset(Japan.data, T1 == 0), 
            "intv.income.gap", "T2", "income.rank.undest", 
            TRUE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by estimate\n(rank)", 
            c("Accurate", "Underest."))
effect.plot(subset(Japan.data, T1 == 0 & income.gap.undest == 1), 
            "intv.income.gap", "T2", "high.income", 
            TRUE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by income\n(Only underest. for ratio)", 
            c("Low", "High"))
effect.plot(subset(Japan.data, T1 == 0 & income.rank.undest == 1), 
            "intv.income.gap", "T2", "high.income", 
            TRUE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by income\n(Only underest. for rank)", 
            c("Low", "High"))
title("Income gap intervention (Japan)", 
      cex.main = 1.5, line = 1, outer = TRUE)
dev.off()

# Figure 4
cairo_pdf("Figure_4_1.pdf", width = 6, height = 1.5, pointsize = 7)
par(mar = c(1, 2.5, 2, 0.5), lwd = 0.5, oma = c(0, 0, 2, 0))
layout(matrix(1:5, 1, 5, byrow = TRUE))
effect.plot(Korea.data, "intv.gender.gap", "T1", 
            FE = FALSE, ylim = c(-0.5, 1), yaxis.labs = seq(-0.5, 1, 0.25), 
            title = "Entire sample")
effect.plot(subset(Korea.data, T2 == 0), 
            "intv.gender.gap", "T1", "gender.gap.undest", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by estimate\n(ratio)", 
            c("Overest.", "Underest."))
effect.plot(subset(Korea.data, T2 == 0), 
            "intv.gender.gap", "T1", "gender.rank.undest", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by estimate\n(rank)", 
            c("Accurate", "Underest."))
effect.plot(subset(Korea.data, T2 == 0 & gender.gap.undest == 1), 
            "intv.gender.gap", "T1", "male", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by gender\n(Only underest. for ratio)", 
            c("Women", "Men"))
effect.plot(subset(Korea.data, T2 == 0 & gender.rank.undest == 1), 
            "intv.gender.gap", "T1", "male", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by gender\n(Only underest. for rank)", 
            c("Women", "Men"))
title("Gender pay gap intervention (Korea)", 
      cex.main = 1.5, line = 1, outer = TRUE)
dev.off()

cairo_pdf("Figure_4_2.pdf", width = 6, height = 1.5, pointsize = 7)
par(mar = c(1, 2.5, 2, 0.5), lwd = 0.5, oma = c(0, 0, 2, 0))
layout(matrix(1:5, 1, 5, byrow = TRUE))
effect.plot(Japan.data, "intv.gender.gap", "T1", 
            FE = TRUE, ylim = c(-0.5, 1), yaxis.labs = seq(-0.5, 1, 0.25), 
            title = "Entire sample")
effect.plot(subset(Japan.data, T2 == 0), 
            "intv.gender.gap", "T1", "gender.gap.undest", 
            TRUE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by estimate\n(ratio)", 
            c("Overest.", "Underest."))
effect.plot(subset(Japan.data, T2 == 0), 
            "intv.gender.gap", "T1", "gender.rank.undest", 
            TRUE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by estimate\n(rank)", 
            c("Accurate", "Underest."))
effect.plot(subset(Japan.data, T2 == 0 & gender.gap.undest == 1), 
            "intv.gender.gap", "T1", "male", 
            TRUE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by gender\n(Only underest. for ratio)", 
            c("Women", "Men"))
effect.plot(subset(Japan.data, T2 == 0 & gender.rank.undest == 1), 
            "intv.gender.gap", "T1", "male", 
            TRUE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by gender\n(Only underest. for rank)", 
            c("Women", "Men"))
title("Gender pay gap intervention (Japan)", 
      cex.main = 1.5, line = 1, outer = TRUE)
dev.off()

#### comparison of corrective effects across issue domains ####
## point estimates
# South Korea
H3.income.gap.intv.gap.Korea <- lm_robust(intv.income.gap ~ (T1 + T2) * income.gap.overest, 
                                          Korea.data)
H3.gender.gap.intv.gap.Korea <- lm_robust(intv.gender.gap ~ (T1 + T2) * gender.gap.overest, 
                                          Korea.data)
round(H3.income.gap.intv.gap.Korea$coefficients["T2"] - 
        H3.gender.gap.intv.gap.Korea$coefficients["T1"], 3)

# Japan
H3.income.gap.intv.gap.Japan <- lm_robust(intv.income.gap ~ (T1 + T2) * income.gap.overest, 
                                          Japan.data)
H3.gender.gap.intv.gap.Japan <- lm_robust(intv.gender.gap ~ (T1 + T2) * gender.gap.overest, 
                                          Japan.data)
round(H3.income.gap.intv.gap.Japan$coefficients["T2"] - 
        H3.gender.gap.intv.gap.Japan$coefficients["T1"], 3)

## bootstrap confidence intervals
# Note: Results for outcomes other than support for policies to reduce 
# the income gap and the gender pay gap are reported in Table A5 in the appendix.

# South Korea
cl <- makeCluster(10)
registerDoParallel(cl)
H3.boot.Korea <- foreach(i = 1:1000, .packages = "estimatr", .combine = "rbind") %dopar% {
  out <- rep(NA, 8)
  set.seed(i)
  boot.data <- Korea.data[sample(1:nrow(Korea.data), nrow(Korea.data), replace = TRUE), ]
  boot.income.gap.intv.gap <- lm_robust(intv.income.gap ~ (T1 + T2) * income.gap.overest, 
                                        boot.data)
  boot.income.rank.intv.gap <- lm_robust(intv.income.gap ~ (T1 + T2) * income.rank.overest, 
                                         boot.data)
  boot.income.gap.intv.policy <- lm_robust(intv.income.policy ~ (T1 + T2) * income.gap.overest, 
                                           boot.data)
  boot.income.rank.intv.policy <- lm_robust(intv.income.policy ~ (T1 + T2) * income.rank.overest, 
                                            boot.data)
  boot.income.gap.feeling.abs <- lm_robust(feeling.poor ~ (T1 + T2) * income.gap.overest, 
                                           boot.data)
  boot.income.rank.feeling.abs <- lm_robust(feeling.poor ~ (T1 + T2) * income.rank.overest, 
                                            boot.data)
  boot.income.gap.feeling.rel <- lm_robust(feeling.income.diff ~ (T1 + T2) * income.gap.overest, 
                                           boot.data)
  boot.income.rank.feeling.rel <- lm_robust(feeling.income.diff ~ (T1 + T2) * income.rank.overest, 
                                            boot.data)
  boot.gender.gap.intv.gap <- lm_robust(intv.gender.gap ~ (T1 + T2) * gender.gap.overest, 
                                        boot.data)
  boot.gender.rank.intv.gap <- lm_robust(intv.gender.gap ~ (T1 + T2) * gender.rank.overest, 
                                         boot.data)
  boot.gender.gap.intv.policy <- lm_robust(intv.gender.policy ~ (T1 + T2) * gender.gap.overest, 
                                           boot.data)
  boot.gender.rank.intv.policy <- lm_robust(intv.gender.policy ~ (T1 + T2) * gender.rank.overest, 
                                            boot.data)
  boot.gender.gap.feeling.abs <- lm_robust(feeling.women ~ (T1 + T2) * gender.gap.overest, 
                                           boot.data)
  boot.gender.rank.feeling.abs <- lm_robust(feeling.women ~ (T1 + T2) * gender.rank.overest, 
                                            boot.data)
  boot.gender.gap.feeling.rel <- lm_robust(feeling.gender.diff ~ (T1 + T2) * gender.gap.overest, 
                                           boot.data)
  boot.gender.rank.feeling.rel <- lm_robust(feeling.gender.diff ~ (T1 + T2) * gender.rank.overest, 
                                            boot.data)
  out[1] <- boot.income.gap.intv.gap$coefficients["T2"] - 
    boot.gender.gap.intv.gap$coefficients["T1"]
  out[2] <- boot.income.rank.intv.gap$coefficients["T2"] - 
    boot.gender.rank.intv.gap$coefficients["T1"]
  out[3] <- boot.income.gap.intv.policy$coefficients["T2"] - 
    boot.gender.gap.intv.policy$coefficients["T1"]
  out[4] <- boot.income.rank.intv.policy$coefficients["T2"] - 
    boot.gender.rank.intv.policy$coefficients["T1"]
  out[5] <- boot.income.gap.feeling.abs$coefficients["T2"] - 
    boot.gender.gap.feeling.abs$coefficients["T1"]
  out[6] <- boot.income.rank.feeling.abs$coefficients["T2"] - 
    boot.gender.rank.feeling.abs$coefficients["T1"]
  out[7] <- boot.income.gap.feeling.rel$coefficients["T2"] - 
    boot.gender.gap.feeling.rel$coefficients["T1"]
  out[8] <- boot.income.rank.feeling.rel$coefficients["T2"] - 
    boot.gender.rank.feeling.rel$coefficients["T1"]
  out
}
stopCluster(cl)

round(apply(H3.boot.Korea, 2, quantile, probs = c(0.025, 0.975)), 3)

# Japan
cl <- makeCluster(10)
registerDoParallel(cl)
H3.boot.Japan <- foreach(i = 1:1000, .packages = "estimatr", .combine = "rbind") %dopar% {
  out <- rep(NA, 8)
  set.seed(i)
  boot.data <- Japan.data[sample(1:nrow(Japan.data), nrow(Japan.data), replace = TRUE), ]
  boot.income.gap.intv.gap <- lm_robust(intv.income.gap ~ (T1 + T2) * income.gap.overest, 
                                        boot.data, fixed_effects = block)
  boot.income.rank.intv.gap <- lm_robust(intv.income.gap ~ (T1 + T2) * income.rank.overest, 
                                         boot.data, fixed_effects = block)
  boot.income.gap.intv.policy <- lm_robust(intv.income.policy ~ (T1 + T2) * income.gap.overest, 
                                           boot.data, fixed_effects = block)
  boot.income.rank.intv.policy <- lm_robust(intv.income.policy ~ (T1 + T2) * income.rank.overest, 
                                            boot.data, fixed_effects = block)
  boot.income.gap.feeling.abs <- lm_robust(feeling.poor ~ (T1 + T2) * income.gap.overest, 
                                           boot.data, fixed_effects = block)
  boot.income.rank.feeling.abs <- lm_robust(feeling.poor ~ (T1 + T2) * income.rank.overest, 
                                            boot.data, fixed_effects = block)
  boot.income.gap.feeling.rel <- lm_robust(feeling.income.diff ~ (T1 + T2) * income.gap.overest, 
                                           boot.data, fixed_effects = block)
  boot.income.rank.feeling.rel <- lm_robust(feeling.income.diff ~ (T1 + T2) * income.rank.overest, 
                                            boot.data, fixed_effects = block)
  boot.gender.gap.intv.gap <- lm_robust(intv.gender.gap ~ (T1 + T2) * gender.gap.overest, 
                                        boot.data, fixed_effects = block)
  boot.gender.rank.intv.gap <- lm_robust(intv.gender.gap ~ (T1 + T2) * gender.rank.overest, 
                                         boot.data, fixed_effects = block)
  boot.gender.gap.intv.policy <- lm_robust(intv.gender.policy ~ (T1 + T2) * gender.gap.overest, 
                                           boot.data, fixed_effects = block)
  boot.gender.rank.intv.policy <- lm_robust(intv.gender.policy ~ (T1 + T2) * gender.rank.overest, 
                                            boot.data, fixed_effects = block)
  boot.gender.gap.feeling.abs <- lm_robust(feeling.women ~ (T1 + T2) * gender.gap.overest, 
                                           boot.data, fixed_effects = block)
  boot.gender.rank.feeling.abs <- lm_robust(feeling.women ~ (T1 + T2) * gender.rank.overest, 
                                            boot.data, fixed_effects = block)
  boot.gender.gap.feeling.rel <- lm_robust(feeling.gender.diff ~ (T1 + T2) * gender.gap.overest, 
                                           boot.data, fixed_effects = block)
  boot.gender.rank.feeling.rel <- lm_robust(feeling.gender.diff ~ (T1 + T2) * gender.rank.overest, 
                                            boot.data, fixed_effects = block)
  out[1] <- boot.income.gap.intv.gap$coefficients["T2"] - 
    boot.gender.gap.intv.gap$coefficients["T1"]
  out[2] <- boot.income.rank.intv.gap$coefficients["T2"] - 
    boot.gender.rank.intv.gap$coefficients["T1"]
  out[3] <- boot.income.gap.intv.policy$coefficients["T2"] - 
    boot.gender.gap.intv.policy$coefficients["T1"]
  out[4] <- boot.income.rank.intv.policy$coefficients["T2"] - 
    boot.gender.rank.intv.policy$coefficients["T1"]
  out[5] <- boot.income.gap.feeling.abs$coefficients["T2"] - 
    boot.gender.gap.feeling.abs$coefficients["T1"]
  out[6] <- boot.income.rank.feeling.abs$coefficients["T2"] - 
    boot.gender.rank.feeling.abs$coefficients["T1"]
  out[7] <- boot.income.gap.feeling.rel$coefficients["T2"] - 
    boot.gender.gap.feeling.rel$coefficients["T1"]
  out[8] <- boot.income.rank.feeling.rel$coefficients["T2"] - 
    boot.gender.rank.feeling.rel$coefficients["T1"]
  out
}
stopCluster(cl)

round(apply(H3.boot.Japan, 2, quantile, probs = c(0.025, 0.975)), 3)

#### cross-national difference in the effect of gender pay gap information on the primary outcome ####
# Figure 5
cairo_pdf("Figure_5.pdf", width = 4, height = 2, pointsize = 7)
par(mar = c(1, 2.5, 2, 0.5), lwd = 0.5)
layout(matrix(1:2, 1, 2))
effect.plot(subset(combined.data, T2 == 0 & gender.gap.undest == 1), 
            "intv.gender.gap", "T1", "Japan", 
            FALSE, c(-0.5, 0.5), seq(-0.5, 0.5, 0.25), 
            "Gender pay gap intervention\n(Only underest. for ratio)", 
            c("Korea", "Japan"), cex.1 = 1.1, cex.2 = 1)
effect.plot(subset(combined.data, T2 == 0 & gender.rank.undest == 1), 
            "intv.gender.gap", "T1", "Japan", 
            FALSE, c(-0.5, 0.5), seq(-0.5, 0.5, 0.25), 
            "Gender pay gap intervention\n(Only underest. for rank)", 
            c("Korea", "Japan"), cex.1 = 1.1, cex.2 = 1)
dev.off()

#### t-tests to compare overestimators and underestimators ####
## South Korea
# Table A1 (A)
t.test.fun(male ~ gender.gap.undest, subset(Korea.data, T1 == 1))
t.test.fun(age ~ gender.gap.undest, subset(Korea.data, T1 == 1))
t.test.fun(educ.cat ~ gender.gap.undest, subset(Korea.data, T1 == 1))
t.test.fun(income ~ gender.gap.undest, subset(Korea.data, T1 == 1))
t.test.fun(ideology ~ gender.gap.undest, subset(Korea.data, T1 == 1))
t.test.fun(HS ~ gender.gap.undest, subset(Korea.data, T1 == 1))
t.test.fun(BS ~ gender.gap.undest, subset(Korea.data, T1 == 1))
t.test.fun(MS ~ gender.gap.undest, subset(Korea.data, T1 == 1))

t.test.fun(male ~ gender.rank.undest, subset(Korea.data, T1 == 1))
t.test.fun(age ~ gender.rank.undest, subset(Korea.data, T1 == 1))
t.test.fun(educ.cat ~ gender.rank.undest, subset(Korea.data, T1 == 1))
t.test.fun(income ~ gender.rank.undest, subset(Korea.data, T1 == 1))
t.test.fun(ideology ~ gender.rank.undest, subset(Korea.data, T1 == 1))
t.test.fun(HS ~ gender.rank.undest, subset(Korea.data, T1 == 1))
t.test.fun(BS ~ gender.rank.undest, subset(Korea.data, T1 == 1))
t.test.fun(MS ~ gender.rank.undest, subset(Korea.data, T1 == 1))

## Japan
# Table A1 (B)
t.test.fun(male ~ gender.gap.undest, subset(Japan.data, T1 == 1))
t.test.fun(age ~ gender.gap.undest, subset(Japan.data, T1 == 1))
t.test.fun(educ.cat ~ gender.gap.undest, subset(Japan.data, T1 == 1))
t.test.fun(income ~ gender.gap.undest, subset(Japan.data, T1 == 1))
t.test.fun(ideology ~ gender.gap.undest, subset(Japan.data, T1 == 1))
t.test.fun(HS ~ gender.gap.undest, subset(Japan.data, T1 == 1))
t.test.fun(BS ~ gender.gap.undest, subset(Japan.data, T1 == 1))
t.test.fun(MS ~ gender.gap.undest, subset(Japan.data, T1 == 1))

t.test.fun(male ~ gender.rank.undest, subset(Japan.data, T1 == 1))
t.test.fun(age ~ gender.rank.undest, subset(Japan.data, T1 == 1))
t.test.fun(educ.cat ~ gender.rank.undest, subset(Japan.data, T1 == 1))
t.test.fun(income ~ gender.rank.undest, subset(Japan.data, T1 == 1))
t.test.fun(ideology ~ gender.rank.undest, subset(Japan.data, T1 == 1))
t.test.fun(HS ~ gender.rank.undest, subset(Japan.data, T1 == 1))
t.test.fun(BS ~ gender.rank.undest, subset(Japan.data, T1 == 1))
t.test.fun(MS ~ gender.rank.undest, subset(Japan.data, T1 == 1))

#### two-way ANOVA for underestimator-by-gender interaction ####
## South Korea
# Table A2 (A)
two.way.anova.fun("age", "gender.gap.undest", "male", subset(Korea.data, T1 == 1))
two.way.anova.fun("educ.cat", "gender.gap.undest", "male", subset(Korea.data, T1 == 1))
two.way.anova.fun("income", "gender.gap.undest", "male", subset(Korea.data, T1 == 1))
two.way.anova.fun("ideology", "gender.gap.undest", "male", subset(Korea.data, T1 == 1))
two.way.anova.fun("HS", "gender.gap.undest", "male", subset(Korea.data, T1 == 1))
two.way.anova.fun("BS", "gender.gap.undest", "male", subset(Korea.data, T1 == 1))
two.way.anova.fun("MS", "gender.gap.undest", "male", subset(Korea.data, T1 == 1))

# Table A2 (B)
two.way.anova.fun("age", "gender.rank.undest", "male", subset(Korea.data, T1 == 1))
two.way.anova.fun("educ.cat", "gender.rank.undest", "male", subset(Korea.data, T1 == 1))
two.way.anova.fun("income", "gender.rank.undest", "male", subset(Korea.data, T1 == 1))
two.way.anova.fun("ideology", "gender.rank.undest", "male", subset(Korea.data, T1 == 1))
two.way.anova.fun("HS", "gender.rank.undest", "male", subset(Korea.data, T1 == 1))
two.way.anova.fun("BS", "gender.rank.undest", "male", subset(Korea.data, T1 == 1))
two.way.anova.fun("MS", "gender.rank.undest", "male", subset(Korea.data, T1 == 1))

## Japan
# Table A3 (A)
two.way.anova.fun("age", "gender.gap.undest", "male", subset(Japan.data, T1 == 1))
two.way.anova.fun("educ.cat", "gender.gap.undest", "male", subset(Japan.data, T1 == 1))
two.way.anova.fun("income", "gender.gap.undest", "male", subset(Japan.data, T1 == 1))
two.way.anova.fun("ideology", "gender.gap.undest", "male", subset(Japan.data, T1 == 1))
two.way.anova.fun("HS", "gender.gap.undest", "male", subset(Japan.data, T1 == 1))
two.way.anova.fun("BS", "gender.gap.undest", "male", subset(Japan.data, T1 == 1))
two.way.anova.fun("MS", "gender.gap.undest", "male", subset(Japan.data, T1 == 1))

# Table A3 (B)
two.way.anova.fun("age", "gender.rank.undest", "male", subset(Japan.data, T1 == 1))
two.way.anova.fun("educ.cat", "gender.rank.undest", "male", subset(Japan.data, T1 == 1))
two.way.anova.fun("income", "gender.rank.undest", "male", subset(Japan.data, T1 == 1))
two.way.anova.fun("ideology", "gender.rank.undest", "male", subset(Japan.data, T1 == 1))
two.way.anova.fun("HS", "gender.rank.undest", "male", subset(Japan.data, T1 == 1))
two.way.anova.fun("BS", "gender.rank.undest", "male", subset(Japan.data, T1 == 1))
two.way.anova.fun("MS", "gender.rank.undest", "male", subset(Japan.data, T1 == 1))

#### covariate balance checks ####
# Table A4
equiv.test.Korea <- 
  rbind(equivalence.F.test(Korea.data, "male", "condition", 0.25), 
        equivalence.F.test(Korea.data, "age", "condition", 0.25), 
        equivalence.F.test(Korea.data, "middle.educ", "condition", 0.25), 
        equivalence.F.test(Korea.data, "high.educ", "condition", 0.25), 
        equivalence.F.test(Korea.data, "income", "condition", 0.25), 
        equivalence.F.test(Korea.data, "ideology", "condition", 0.25), 
        equivalence.F.test(Korea.data, "HS", "condition", 0.25), 
        equivalence.F.test(Korea.data, "BS", "condition", 0.25), 
        equivalence.F.test(Korea.data, "MS", "condition", 0.25), 
        equivalence.F.test(Korea.data, "income.gap.est", "condition", 0.25), 
        equivalence.F.test(Korea.data, "income.rank.est", "condition", 0.25), 
        equivalence.F.test(Korea.data, "gender.gap.est", "condition", 0.25), 
        equivalence.F.test(Korea.data, "gender.rank.est", "condition", 0.25))
equiv.test.Korea

equiv.test.Japan <- 
  rbind(equivalence.F.test(Japan.data, "male", "condition", 0.25), 
        equivalence.F.test(Japan.data, "age", "condition", 0.25), 
        equivalence.F.test(Japan.data, "middle.educ", "condition", 0.25), 
        equivalence.F.test(Japan.data, "high.educ", "condition", 0.25), 
        equivalence.F.test(Japan.data, "income", "condition", 0.25), 
        equivalence.F.test(Japan.data, "ideology", "condition", 0.25), 
        equivalence.F.test(Japan.data, "HS", "condition", 0.25), 
        equivalence.F.test(Japan.data, "BS", "condition", 0.25), 
        equivalence.F.test(Japan.data, "MS", "condition", 0.25), 
        equivalence.F.test(Japan.data, "income.gap.est", "condition", 0.25), 
        equivalence.F.test(Japan.data, "income.rank.est", "condition", 0.25), 
        equivalence.F.test(Japan.data, "gender.gap.est", "condition", 0.25), 
        equivalence.F.test(Japan.data, "gender.rank.est", "condition", 0.25))
equiv.test.Japan

#### results for the other outcomes ####
# Figure A1
cairo_pdf("Figure_A1_1.pdf", width = 6, height = 1.5, pointsize = 7)
par(mar = c(1, 2.5, 2, 0.5), lwd = 0.5, oma = c(0, 0, 2, 0))
layout(matrix(1:5, 1, 5, byrow = TRUE))
effect.plot(Korea.data, "intv.income.policy", "T2", 
            FE = FALSE, ylim = c(-0.5, 1), yaxis.labs = seq(-0.5, 1, 0.25), 
            title = "Entire sample")
effect.plot(subset(Korea.data, T1 == 0), 
            "intv.income.policy", "T2", "income.gap.undest", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by estimate\n(ratio)", 
            c("Overest.", "Underest."))
effect.plot(subset(Korea.data, T1 == 0), 
            "intv.income.policy", "T2", "income.rank.undest", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by estimate\n(rank)", 
            c("Accurate", "Underest."))
effect.plot(subset(Korea.data, T1 == 0 & income.gap.undest == 1), 
            "intv.income.policy", "T2", "high.income", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by income\n(Only underest. for ratio)", 
            c("Low", "High"))
effect.plot(subset(Korea.data, T1 == 0 & income.rank.undest == 1), 
            "intv.income.policy", "T2", "high.income", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by income\n(Only underest. for rank)", 
            c("Low", "High"))
title("Welfare for poor people (Korea)", 
      cex.main = 1.5, line = 1, outer = TRUE)
dev.off()

cairo_pdf("Figure_A1_2.pdf", width = 6, height = 1.5, pointsize = 7)
par(mar = c(1, 2.5, 2, 0.5), lwd = 0.5, oma = c(0, 0, 2, 0))
layout(matrix(1:5, 1, 5, byrow = TRUE))
effect.plot(Japan.data, "intv.income.policy", "T2", 
            FE = TRUE, ylim = c(-0.5, 1), yaxis.labs = seq(-0.5, 1, 0.25), 
            title = "Entire sample")
effect.plot(subset(Japan.data, T1 == 0), 
            "intv.income.policy", "T2", "income.gap.undest", 
            TRUE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by estimate\n(ratio)", 
            c("Overest.", "Underest."))
effect.plot(subset(Japan.data, T1 == 0), 
            "intv.income.policy", "T2", "income.rank.undest", 
            TRUE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by estimate\n(rank)", 
            c("Accurate", "Underest."))
effect.plot(subset(Japan.data, T1 == 0 & income.gap.undest == 1), 
            "intv.income.policy", "T2", "high.income", 
            TRUE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by income\n(Only underest. for ratio)", 
            c("Low", "High"))
effect.plot(subset(Japan.data, T1 == 0 & income.rank.undest == 1), 
            "intv.income.policy", "T2", "high.income", 
            TRUE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by income\n(Only underest. for rank)", 
            c("Low", "High"))
title("Welfare for poor people (Japan)", 
      cex.main = 1.5, line = 1, outer = TRUE)
dev.off()

# Figure A2
cairo_pdf("Figure_A2_1.pdf", width = 6, height = 1.5, pointsize = 7)
par(mar = c(1, 2.5, 2, 0.5), lwd = 0.5, oma = c(0, 0, 2, 0))
layout(matrix(1:5, 1, 5, byrow = TRUE))
effect.plot(Korea.data, "intv.gender.policy", "T1", 
            FE = FALSE, ylim = c(-0.5, 1), yaxis.labs = seq(-0.5, 1, 0.25), 
            title = "Entire sample")
effect.plot(subset(Korea.data, T2 == 0), 
            "intv.gender.policy", "T1", "gender.gap.undest", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by estimate\n(ratio)", 
            c("Overest.", "Underest."))
effect.plot(subset(Korea.data, T2 == 0), 
            "intv.gender.policy", "T1", "gender.rank.undest", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by estimate\n(rank)", 
            c("Accurate", "Underest."))
effect.plot(subset(Korea.data, T2 == 0 & gender.gap.undest == 1), 
            "intv.gender.policy", "T1", "male", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by gender\n(Only underest. for ratio)", 
            c("Women", "Men"))
effect.plot(subset(Korea.data, T2 == 0 & gender.rank.undest == 1), 
            "intv.gender.policy", "T1", "male", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by gender\n(Only underest. for rank)", 
            c("Women", "Men"))
title("Equal Employment Act (Korea)", 
      cex.main = 1.5, line = 1, outer = TRUE)
dev.off()

cairo_pdf("Figure_A2_2.pdf", width = 6, height = 1.5, pointsize = 7)
par(mar = c(1, 2.5, 2, 0.5), lwd = 0.5, oma = c(0, 0, 2, 0))
layout(matrix(1:5, 1, 5, byrow = TRUE))
effect.plot(Japan.data, "intv.gender.policy", "T1", 
            FE = TRUE, ylim = c(-0.5, 1), yaxis.labs = seq(-0.5, 1, 0.25), 
            title = "Entire sample")
effect.plot(subset(Japan.data, T2 == 0), 
            "intv.gender.policy", "T1", "gender.gap.undest", 
            TRUE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by estimate\n(ratio)", 
            c("Overest.", "Underest."))
effect.plot(subset(Japan.data, T2 == 0), 
            "intv.gender.policy", "T1", "gender.rank.undest", 
            TRUE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by estimate\n(rank)", 
            c("Accurate", "Underest."))
effect.plot(subset(Japan.data, T2 == 0 & gender.gap.undest == 1), 
            "intv.gender.policy", "T1", "male", 
            TRUE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by gender\n(Only underest. for ratio)", 
            c("Women", "Men"))
effect.plot(subset(Japan.data, T2 == 0 & gender.rank.undest == 1), 
            "intv.gender.policy", "T1", "male", 
            TRUE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Heterogeneity by gender\n(Only underest. for rank)", 
            c("Women", "Men"))
title("Equal Employment Act (Japan)", 
      cex.main = 1.5, line = 1, outer = TRUE)
dev.off()

# Figure A3
cairo_pdf("Figure_A3_1.pdf", width = 6, height = 1.5, pointsize = 7)
par(mar = c(1, 2.5, 2, 0.5), lwd = 0.5, oma = c(0, 0, 2, 0))
layout(matrix(1:5, 1, 5, byrow = TRUE))
effect.plot(Korea.data, "feeling.poor", "T2", 
            FE = FALSE, ylim = c(-10, 10), yaxis.labs = seq(-10, 10, 5), 
            title = "Entire sample")
effect.plot(subset(Korea.data, T1 == 0), 
            "feeling.poor", "T2", "income.gap.undest", 
            FALSE, c(-10, 10), seq(-10, 10, 5), 
            "Heterogeneity by estimate\n(ratio)", 
            c("Overest.", "Underest."))
effect.plot(subset(Korea.data, T1 == 0), 
            "feeling.poor", "T2", "income.rank.undest", 
            FALSE, c(-10, 10), seq(-10, 10, 5), 
            "Heterogeneity by estimate\n(rank)", 
            c("Accurate", "Underest."))
effect.plot(subset(Korea.data, T1 == 0 & income.gap.undest == 1), 
            "feeling.poor", "T2", "male", 
            FALSE, c(-10, 10), seq(-10, 10, 5), 
            "Heterogeneity by income\n(Only underest. for ratio)", 
            c("Low", "High"))
effect.plot(subset(Korea.data, T1 == 0 & income.rank.undest == 1), 
            "feeling.poor", "T2", "high.income", 
            FALSE, c(-10, 10), seq(-10, 10, 5), 
            "Heterogeneity by income\n(Only underest. for rank)", 
            c("Low", "High"))
title("Feeling thermometer toward poor people (Korea)", 
      cex.main = 1.5, line = 1, outer = TRUE)
dev.off()

cairo_pdf("Figure_A3_2.pdf", width = 6, height = 1.5, pointsize = 7)
par(mar = c(1, 2.5, 2, 0.5), lwd = 0.5, oma = c(0, 0, 2, 0))
layout(matrix(1:5, 1, 5, byrow = TRUE))
effect.plot(Japan.data, "feeling.poor", "T2", 
            FE = TRUE, ylim = c(-10, 10), yaxis.labs = seq(-10, 10, 5), 
            title = "Entire sample")
effect.plot(subset(Japan.data, T1 == 0), 
            "feeling.poor", "T2", "income.gap.undest", 
            TRUE, c(-10, 10), seq(-10, 10, 5), 
            "Heterogeneity by estimate\n(ratio)", 
            c("Overest.", "Underest."))
effect.plot(subset(Japan.data, T1 == 0), 
            "feeling.poor", "T2", "income.rank.undest", 
            TRUE, c(-10, 10), seq(-10, 10, 5), 
            "Heterogeneity by estimate\n(rank)", 
            c("Accurate", "Underest."))
effect.plot(subset(Japan.data, T1 == 0 & income.gap.undest == 1), 
            "feeling.poor", "T2", "male", 
            TRUE, c(-10, 10), seq(-10, 10, 5), 
            "Heterogeneity by income\n(Only underest. for ratio)", 
            c("Low", "High"))
effect.plot(subset(Japan.data, T1 == 0 & income.rank.undest == 1), 
            "feeling.poor", "T2", "high.income", 
            TRUE, c(-10, 10), seq(-10, 10, 5), 
            "Heterogeneity by income\n(Only underest. for rank)", 
            c("Low", "High"))
title("Feeling thermometer toward poor people (Japan)", 
      cex.main = 1.5, line = 1, outer = TRUE)
dev.off()

# Figure A4
cairo_pdf("Figure_A4_1.pdf", width = 6, height = 1.5, pointsize = 7)
par(mar = c(1, 2.5, 2, 0.5), lwd = 0.5, oma = c(0, 0, 2, 0))
layout(matrix(1:5, 1, 5, byrow = TRUE))
effect.plot(Korea.data, "feeling.women", "T1", 
            FE = FALSE, ylim = c(-10, 10), yaxis.labs = seq(-10, 10, 5), 
            title = "Entire sample")
effect.plot(subset(Korea.data, T2 == 0), 
            "feeling.women", "T1", "gender.gap.undest", 
            FALSE, c(-10, 10), seq(-10, 10, 5), 
            "Heterogeneity by estimate\n(ratio)", 
            c("Overest.", "Underest."))
effect.plot(subset(Korea.data, T2 == 0), 
            "feeling.women", "T1", "gender.rank.undest", 
            FALSE, c(-10, 10), seq(-10, 10, 5), 
            "Heterogeneity by estimate\n(rank)", 
            c("Accurate", "Underest."))
effect.plot(subset(Korea.data, T2 == 0 & gender.gap.undest == 1), 
            "feeling.women", "T1", "male", 
            FALSE, c(-10, 10), seq(-10, 10, 5), 
            "Heterogeneity by gender\n(Only underest. for ratio)", 
            c("Women", "Men"))
effect.plot(subset(Korea.data, T2 == 0 & gender.rank.undest == 1), 
            "feeling.women", "T1", "male", 
            FALSE, c(-10, 10), seq(-10, 10, 5), 
            "Heterogeneity by gender\n(Only underest. for rank)", 
            c("Women", "Men"))
title("Feeling thermometer toward women (Korea)", 
      cex.main = 1.5, line = 1, outer = TRUE)
dev.off()

cairo_pdf("Figure_A4_2.pdf", width = 6, height = 1.5, pointsize = 7)
par(mar = c(1, 2.5, 2, 0.5), lwd = 0.5, oma = c(0, 0, 2, 0))
layout(matrix(1:5, 1, 5, byrow = TRUE))
effect.plot(Japan.data, "feeling.women", "T1", 
            FE = TRUE, ylim = c(-10, 10), yaxis.labs = seq(-10, 10, 5), 
            title = "Entire sample")
effect.plot(subset(Japan.data, T2 == 0), 
            "feeling.women", "T1", "gender.gap.undest", 
            TRUE, c(-10, 10), seq(-10, 10, 5), 
            "Heterogeneity by estimate\n(ratio)", 
            c("Overest.", "Underest."))
effect.plot(subset(Japan.data, T2 == 0), 
            "feeling.women", "T1", "gender.rank.undest", 
            TRUE, c(-10, 10), seq(-10, 10, 5), 
            "Heterogeneity by estimate\n(rank)", 
            c("Accurate", "Underest."))
effect.plot(subset(Japan.data, T2 == 0 & gender.gap.undest == 1), 
            "feeling.women", "T1", "male", 
            TRUE, c(-10, 10), seq(-10, 10, 5), 
            "Heterogeneity by gender\n(Only underest. for ratio)", 
            c("Women", "Men"))
effect.plot(subset(Japan.data, T2 == 0 & gender.rank.undest == 1), 
            "feeling.women", "T1", "male", 
            TRUE, c(-10, 10), seq(-10, 10, 5), 
            "Heterogeneity by gender\n(Only underest. for rank)", 
            c("Women", "Men"))
title("Feeling thermometer toward women (Japan)", 
      cex.main = 1.5, line = 1, outer = TRUE)
dev.off()

# Figure A5
cairo_pdf("Figure_A5_1.pdf", width = 6, height = 1.5, pointsize = 7)
par(mar = c(1, 2.5, 2, 0.5), lwd = 0.5, oma = c(0, 0, 2, 0))
layout(matrix(1:5, 1, 5, byrow = TRUE))
effect.plot(Korea.data, "feeling.income.diff", "T2", 
            FE = FALSE, ylim = c(-12, 12), yaxis.labs = seq(-12, 12, 4), 
            title = "Entire sample")
effect.plot(subset(Korea.data, T1 == 0), 
            "feeling.income.diff", "T2", "income.gap.undest", 
            FALSE, c(-12, 12), seq(-12, 12, 4), 
            "Heterogeneity by estimate\n(ratio)", 
            c("Overest.", "Underest."))
effect.plot(subset(Korea.data, T1 == 0), 
            "feeling.income.diff", "T2", "income.rank.undest", 
            FALSE, c(-12, 12), seq(-12, 12, 4), 
            "Heterogeneity by estimate\n(rank)", 
            c("Accurate", "Underest."))
effect.plot(subset(Korea.data, T1 == 0 & income.gap.undest == 1), 
            "feeling.income.diff", "T2", "male", 
            FALSE, c(-12, 12), seq(-12, 12, 4), 
            "Heterogeneity by income\n(Only underest. for ratio)", 
            c("Low", "High"))
effect.plot(subset(Korea.data, T1 == 0 & income.rank.undest == 1), 
            "feeling.income.diff", "T2", "high.income", 
            FALSE, c(-12, 12), seq(-12, 12, 4), 
            "Heterogeneity by income\n(Only underest. for rank)", 
            c("Low", "High"))
title("Difference between feeling thermometers toward poor and rich (Korea)", 
      cex.main = 1.5, line = 1, outer = TRUE)
dev.off()

cairo_pdf("Figure_A5_2.pdf", width = 6, height = 1.5, pointsize = 7)
par(mar = c(1, 2.5, 2, 0.5), lwd = 0.5, oma = c(0, 0, 2, 0))
layout(matrix(1:5, 1, 5, byrow = TRUE))
effect.plot(Japan.data, "feeling.income.diff", "T2", 
            FE = TRUE, ylim = c(-12, 12), yaxis.labs = seq(-12, 12, 4), 
            title = "Entire sample")
effect.plot(subset(Japan.data, T1 == 0), 
            "feeling.income.diff", "T2", "income.gap.undest", 
            TRUE, c(-12, 12), seq(-12, 12, 4), 
            "Heterogeneity by estimate\n(ratio)", 
            c("Overest.", "Underest."))
effect.plot(subset(Japan.data, T1 == 0), 
            "feeling.income.diff", "T2", "income.rank.undest", 
            TRUE, c(-12, 12), seq(-12, 12, 4), 
            "Heterogeneity by estimate\n(rank)", 
            c("Accurate", "Underest."))
effect.plot(subset(Japan.data, T1 == 0 & income.gap.undest == 1), 
            "feeling.income.diff", "T2", "male", 
            TRUE, c(-12, 12), seq(-12, 12, 4), 
            "Heterogeneity by income\n(Only underest. for ratio)", 
            c("Low", "High"))
effect.plot(subset(Japan.data, T1 == 0 & income.rank.undest == 1), 
            "feeling.income.diff", "T2", "high.income", 
            TRUE, c(-12, 12), seq(-12, 12, 4), 
            "Heterogeneity by income\n(Only underest. for rank)", 
            c("Low", "High"))
title("Difference between feeling thermometers toward poor and rich (Japan)", 
      cex.main = 1.5, line = 1, outer = TRUE)
dev.off()

# Figure A6
cairo_pdf("Figure_A6_1.pdf", width = 6, height = 1.5, pointsize = 7)
par(mar = c(1, 2.5, 2, 0.5), lwd = 0.5, oma = c(0, 0, 2, 0))
layout(matrix(1:5, 1, 5, byrow = TRUE))
effect.plot(Korea.data, "feeling.gender.diff", "T1", 
            FE = FALSE, ylim = c(-12, 12), yaxis.labs = seq(-12, 12, 4), 
            title = "Entire sample")
effect.plot(subset(Korea.data, T2 == 0), 
            "feeling.gender.diff", "T1", "gender.gap.undest", 
            FALSE, c(-12, 12), seq(-12, 12, 4), 
            "Heterogeneity by estimate\n(ratio)", 
            c("Overest.", "Underest."))
effect.plot(subset(Korea.data, T2 == 0), 
            "feeling.gender.diff", "T1", "gender.rank.undest", 
            FALSE, c(-12, 12), seq(-12, 12, 4), 
            "Heterogeneity by estimate\n(rank)", 
            c("Accurate", "Underest."))
effect.plot(subset(Korea.data, T2 == 0 & gender.gap.undest == 1), 
            "feeling.gender.diff", "T1", "male", 
            FALSE, c(-12, 12), seq(-12, 12, 4), 
            "Heterogeneity by gender\n(Only underest. for ratio)", 
            c("Women", "Men"))
effect.plot(subset(Korea.data, T2 == 0 & gender.rank.undest == 1), 
            "feeling.gender.diff", "T1", "male", 
            FALSE, c(-12, 12), seq(-12, 12, 4), 
            "Heterogeneity by gender\n(Only underest. for rank)", 
            c("Women", "Men"))
title("Difference between feeling thermometers toward women and men (Korea)", 
      cex.main = 1.5, line = 1, outer = TRUE)
dev.off()

cairo_pdf("Figure_A6_2.pdf", width = 6, height = 1.5, pointsize = 7)
par(mar = c(1, 2.5, 2, 0.5), lwd = 0.5, oma = c(0, 0, 2, 0))
layout(matrix(1:5, 1, 5, byrow = TRUE))
effect.plot(Japan.data, "feeling.gender.diff", "T1", 
            FE = TRUE, ylim = c(-12, 12), yaxis.labs = seq(-12, 12, 4), 
            title = "Entire sample")
effect.plot(subset(Japan.data, T2 == 0), 
            "feeling.gender.diff", "T1", "gender.gap.undest", 
            TRUE, c(-12, 12), seq(-12, 12, 4), 
            "Heterogeneity by estimate\n(ratio)", 
            c("Overest.", "Underest."))
effect.plot(subset(Japan.data, T2 == 0), 
            "feeling.gender.diff", "T1", "gender.rank.undest", 
            TRUE, c(-12, 12), seq(-12, 12, 4), 
            "Heterogeneity by estimate\n(rank)", 
            c("Accurate", "Underest."))
effect.plot(subset(Japan.data, T2 == 0 & gender.gap.undest == 1), 
            "feeling.gender.diff", "T1", "male", 
            TRUE, c(-12, 12), seq(-12, 12, 4), 
            "Heterogeneity by gender\n(Only underest. for ratio)", 
            c("Women", "Men"))
effect.plot(subset(Japan.data, T2 == 0 & gender.rank.undest == 1), 
            "feeling.gender.diff", "T1", "male", 
            TRUE, c(-12, 12), seq(-12, 12, 4), 
            "Heterogeneity by gender\n(Only underest. for rank)", 
            c("Women", "Men"))
title("Difference between feeling thermometers toward women and men (Japan)", 
      cex.main = 1.5, line = 1, outer = TRUE)
dev.off()

# Table A5
# Note: A point estimate in the left-most column and bootstrap confidence 
# intervals in all columns are produced in the "comparison of corrective 
# effects across issue domains" section
H3.income.rank.intv.gap.Korea <- lm_robust(intv.income.gap ~ (T1 + T2) * income.rank.overest, 
                                          Korea.data)
H3.gender.rank.intv.gap.Korea <- lm_robust(intv.gender.gap ~ (T1 + T2) * gender.rank.overest, 
                                          Korea.data)
round(H3.income.rank.intv.gap.Korea$coefficients["T2"] - 
        H3.gender.rank.intv.gap.Korea$coefficients["T1"], 3)

H3.income.gap.intv.policy.Korea <- lm_robust(intv.income.policy ~ (T1 + T2) * income.gap.overest, 
                                             Korea.data)
H3.gender.gap.intv.policy.Korea <- lm_robust(intv.gender.policy ~ (T1 + T2) * gender.gap.overest, 
                                             Korea.data)
round(H3.income.gap.intv.policy.Korea$coefficients["T2"] - 
        H3.gender.gap.intv.policy.Korea$coefficients["T1"], 3)

H3.income.rank.intv.policy.Korea <- lm_robust(intv.income.policy ~ (T1 + T2) * income.rank.overest, 
                                             Korea.data)
H3.gender.rank.intv.policy.Korea <- lm_robust(intv.gender.policy ~ (T1 + T2) * gender.rank.overest, 
                                             Korea.data)
round(H3.income.rank.intv.policy.Korea$coefficients["T2"] - 
        H3.gender.rank.intv.policy.Korea$coefficients["T1"], 3)

H3.income.gap.feeling.abs.Korea <- lm_robust(feeling.poor ~ (T1 + T2) * income.gap.overest, 
                                             Korea.data)
H3.gender.gap.feeling.abs.Korea <- lm_robust(feeling.women ~ (T1 + T2) * gender.gap.overest, 
                                             Korea.data)
round(H3.income.gap.feeling.abs.Korea$coefficients["T2"] - 
        H3.gender.gap.feeling.abs.Korea$coefficients["T1"], 3)

H3.income.rank.feeling.abs.Korea <- lm_robust(feeling.poor ~ (T1 + T2) * income.rank.overest, 
                                              Korea.data)
H3.gender.rank.feeling.abs.Korea <- lm_robust(feeling.women ~ (T1 + T2) * gender.rank.overest, 
                                              Korea.data)
round(H3.income.rank.feeling.abs.Korea$coefficients["T2"] - 
        H3.gender.rank.feeling.abs.Korea$coefficients["T1"], 3)

H3.income.gap.feeling.rel.Korea <- lm_robust(feeling.income.diff ~ (T1 + T2) * income.gap.overest, 
                                             Korea.data)
H3.gender.gap.feeling.rel.Korea <- lm_robust(feeling.gender.diff ~ (T1 + T2) * gender.gap.overest, 
                                             Korea.data)
round(H3.income.gap.feeling.rel.Korea$coefficients["T2"] - 
        H3.gender.gap.feeling.rel.Korea$coefficients["T1"], 3)

H3.income.rank.feeling.rel.Korea <- lm_robust(feeling.income.diff ~ (T1 + T2) * income.rank.overest, 
                                              Korea.data)
H3.gender.rank.feeling.rel.Korea <- lm_robust(feeling.gender.diff ~ (T1 + T2) * gender.rank.overest, 
                                              Korea.data)
round(H3.income.rank.feeling.rel.Korea$coefficients["T2"] - 
        H3.gender.rank.feeling.rel.Korea$coefficients["T1"], 3)

H3.income.rank.intv.gap.Japan <- lm_robust(intv.income.gap ~ (T1 + T2) * income.rank.overest, 
                                           Japan.data)
H3.gender.rank.intv.gap.Japan <- lm_robust(intv.gender.gap ~ (T1 + T2) * gender.rank.overest, 
                                           Japan.data)
round(H3.income.rank.intv.gap.Japan$coefficients["T2"] - 
        H3.gender.rank.intv.gap.Japan$coefficients["T1"], 3)

H3.income.gap.intv.policy.Japan <- lm_robust(intv.income.policy ~ (T1 + T2) * income.gap.overest, 
                                             Japan.data)
H3.gender.gap.intv.policy.Japan <- lm_robust(intv.gender.policy ~ (T1 + T2) * gender.gap.overest, 
                                             Japan.data)
round(H3.income.gap.intv.policy.Japan$coefficients["T2"] - 
        H3.gender.gap.intv.policy.Japan$coefficients["T1"], 3)

H3.income.rank.intv.policy.Japan <- lm_robust(intv.income.policy ~ (T1 + T2) * income.rank.overest, 
                                              Japan.data)
H3.gender.rank.intv.policy.Japan <- lm_robust(intv.gender.policy ~ (T1 + T2) * gender.rank.overest, 
                                              Japan.data)
round(H3.income.rank.intv.policy.Japan$coefficients["T2"] - 
        H3.gender.rank.intv.policy.Japan$coefficients["T1"], 3)

H3.income.gap.feeling.abs.Japan <- lm_robust(feeling.poor ~ (T1 + T2) * income.gap.overest, 
                                             Japan.data)
H3.gender.gap.feeling.abs.Japan <- lm_robust(feeling.women ~ (T1 + T2) * gender.gap.overest, 
                                             Japan.data)
round(H3.income.gap.feeling.abs.Japan$coefficients["T2"] - 
        H3.gender.gap.feeling.abs.Japan$coefficients["T1"], 3)

H3.income.rank.feeling.abs.Japan <- lm_robust(feeling.poor ~ (T1 + T2) * income.rank.overest, 
                                              Japan.data)
H3.gender.rank.feeling.abs.Japan <- lm_robust(feeling.women ~ (T1 + T2) * gender.rank.overest, 
                                              Japan.data)
round(H3.income.rank.feeling.abs.Japan$coefficients["T2"] - 
        H3.gender.rank.feeling.abs.Japan$coefficients["T1"], 3)

H3.income.gap.feeling.rel.Japan <- lm_robust(feeling.income.diff ~ (T1 + T2) * income.gap.overest, 
                                             Japan.data)
H3.gender.gap.feeling.rel.Japan <- lm_robust(feeling.gender.diff ~ (T1 + T2) * gender.gap.overest, 
                                             Japan.data)
round(H3.income.gap.feeling.rel.Japan$coefficients["T2"] - 
        H3.gender.gap.feeling.rel.Japan$coefficients["T1"], 3)

H3.income.rank.feeling.rel.Japan <- lm_robust(feeling.income.diff ~ (T1 + T2) * income.rank.overest, 
                                              Japan.data)
H3.gender.rank.feeling.rel.Japan <- lm_robust(feeling.gender.diff ~ (T1 + T2) * gender.rank.overest, 
                                              Japan.data)
round(H3.income.rank.feeling.rel.Japan$coefficients["T2"] - 
        H3.gender.rank.feeling.rel.Japan$coefficients["T1"], 3)

# Figure A7
cairo_pdf("Figure_A7_1.pdf", 
          width = 4.8, height = 1.4, pointsize = 7)
par(mar = c(1, 2.5, 4, 0.5), lwd = 0.5)
layout(matrix(1:4, 1, 4, byrow = TRUE))
effect.plot(subset(combined.data, T2 == 0 & gender.gap.undest == 1), 
            "intv.gender.gap", "T1", "Japan", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Gender pay gap intervention\n(Only underest. for ratio)", 
            c("Korea", "Japan"))
effect.plot(subset(combined.data, T2 == 0 & gender.gap.undest == 1), 
            "intv.gender.policy", "T1", "Japan", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Equal Employment Act\n(Only underest. for ratio)", 
            c("Korea", "Japan"))
effect.plot(subset(combined.data, T2 == 0 & gender.gap.undest == 1), 
            "feeling.women", "T1", "Japan", 
            FALSE, c(-10, 10), seq(-10, 10, 5), 
            "Feeling twd. women\n(Only underest. for ratio)", 
            c("Korea", "Japan"))
effect.plot(subset(combined.data, T2 == 0 & gender.gap.undest == 1), 
            "feeling.gender.diff", "T1", "Japan", 
            FALSE, c(-12, 12), seq(-12, 12, 4), 
            "Diff. b/w feeling\n(Only underest. for ratio)", 
            c("Korea", "Japan"))
dev.off()

cairo_pdf("Figure_A7_2.pdf", 
          width = 4.8, height = 1.4, pointsize = 7)
par(mar = c(1, 2.5, 4, 0.5), lwd = 0.5)
layout(matrix(1:4, 1, 4, byrow = TRUE))
effect.plot(subset(combined.data, T2 == 0 & gender.rank.undest == 1), 
            "intv.gender.gap", "T1", "Japan", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Gender pay gap intervention\n(Only underest. for rank)", 
            c("Korea", "Japan"))
effect.plot(subset(combined.data, T2 == 0 & gender.rank.undest == 1), 
            "intv.gender.policy", "T1", "Japan", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Equal Employment Act\n(Only underest. for rank)", 
            c("Korea", "Japan"))
effect.plot(subset(combined.data, T2 == 0 & gender.rank.undest == 1), 
            "feeling.women", "T1", "Japan", 
            FALSE, c(-10, 10), seq(-10, 10, 5), 
            "Feeling twd. women\n(Only underest. for rank)", 
            c("Korea", "Japan"))
effect.plot(subset(combined.data, T2 == 0 & gender.rank.undest == 1), 
            "feeling.gender.diff", "T1", "Japan", 
            FALSE, c(-12, 12), seq(-12, 12, 4), 
            "Diff. b/w feeling\n(Only underest. for rank)", 
            c("Korea", "Japan"))
dev.off()

#### H4-H6 analyses using all respondents ####
# Figure A8
cairo_pdf("Figure_A8_1.pdf", width = 4.8, height = 1.4, pointsize = 7)
par(mar = c(1, 2.5, 4, 0.5), lwd = 0.5, oma = c(0, 0, 1, 0))
layout(matrix(1:4, 1, 4, byrow = TRUE))
effect.plot(subset(Korea.data, T1 == 0), 
            "intv.income.gap", "T2", "high.income", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Income gap intervention\n", 
            c("Low", "High"))
effect.plot(subset(Korea.data, T1 == 0), 
            "intv.income.policy", "T2", "high.income", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Welfare for poor people\n", 
            c("Low", "High"))
effect.plot(subset(Korea.data, T1 == 0), 
            "feeling.poor", "T2", "high.income", 
            FALSE, c(-10, 10), seq(-10, 10, 5), 
            "Feeling thermometer\ntoward poor people", 
            c("Low", "High"))
effect.plot(subset(Korea.data, T1 == 0), 
            "feeling.income.diff", "T2", "high.income", 
            FALSE, c(-12, 12), seq(-12, 12, 4), 
            "Difference between\nfeeling thermometers\ntoward poor and rich", 
            c("Low", "High"))
title("Heterogeneous treatment effect by income, using all respondents (Korea)", 
      cex.main = 1.5, outer = TRUE)
dev.off()

cairo_pdf("Figure_A8_2.pdf", 
          width = 4.8, height = 1.4, pointsize = 7)
par(mar = c(1, 2.5, 4, 0.5), lwd = 0.5, oma = c(0, 0, 1, 0))
layout(matrix(1:4, 1, 4, byrow = TRUE))
effect.plot(subset(Japan.data, T1 == 0), 
            "intv.income.gap", "T2", "high.income", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Income gap intervention\n", 
            c("Low", "High"))
effect.plot(subset(Japan.data, T1 == 0), 
            "intv.income.policy", "T2", "high.income", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Welfare for poor people\n", 
            c("Low", "High"))
effect.plot(subset(Japan.data, T1 == 0), 
            "feeling.poor", "T2", "high.income", 
            FALSE, c(-10, 10), seq(-10, 10, 5), 
            "Feeling thermometer\ntoward poor people", 
            c("Low", "High"))
effect.plot(subset(Japan.data, T1 == 0), 
            "feeling.income.diff", "T2", "high.income", 
            FALSE, c(-12, 12), seq(-12, 12, 4), 
            "Difference between\nfeeling thermometers\ntoward poor and rich", 
            c("Low", "High"))
title("Heterogeneous treatment effect by income, using all respondents (Japan)", 
      cex.main = 1.5, outer = TRUE)
dev.off()

# Figure A9
cairo_pdf("Figure_A9_1.pdf", width = 4.8, height = 1.4, pointsize = 7)
par(mar = c(1, 2.5, 4, 0.5), lwd = 0.5, oma = c(0, 0, 1, 0))
layout(matrix(1:4, 1, 4, byrow = TRUE))
effect.plot(subset(Korea.data, T2 == 0), 
            "intv.gender.gap", "T1", "male", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Gender pay gap\nintervention", 
            c("Women", "Men"))
effect.plot(subset(Korea.data, T2 == 0), 
            "intv.gender.policy", "T1", "male", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Equal Employment Act\n", 
            c("Women", "Men"))
effect.plot(subset(Korea.data, T2 == 0), 
            "feeling.women", "T1", "male", 
            FALSE, c(-10, 10), seq(-10, 10, 5), 
            "Feeling thermometer\ntoward women", 
            c("Women", "Men"))
effect.plot(subset(Korea.data, T2 == 0), 
            "feeling.gender.diff", "T1", "male", 
            FALSE, c(-12, 12), seq(-12, 12, 4), 
            "Difference between\nfeeling thermometers\ntoward women and men", 
            c("Women", "Men"))
title("Heterogeneous treatment effect by gender, using all respondents (Korea)", 
      cex.main = 1.5, outer = TRUE)
dev.off()

cairo_pdf("Figure_A9_2.pdf", width = 4.8, height = 1.4, pointsize = 7)
par(mar = c(1, 2.5, 4, 0.5), lwd = 0.5, oma = c(0, 0, 1, 0))
layout(matrix(1:4, 1, 4, byrow = TRUE))
effect.plot(subset(Japan.data, T2 == 0), 
            "intv.gender.gap", "T1", "male", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Gender pay gap\nintervention", 
            c("Women", "Men"))
effect.plot(subset(Japan.data, T2 == 0), 
            "intv.gender.policy", "T1", "male", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Equal Employment Act\n", 
            c("Women", "Men"))
effect.plot(subset(Japan.data, T2 == 0), 
            "feeling.women", "T1", "male", 
            FALSE, c(-10, 10), seq(-10, 10, 5), 
            "Feeling thermometer\ntoward women", 
            c("Women", "Men"))
effect.plot(subset(Japan.data, T2 == 0), 
            "feeling.gender.diff", "T1", "male", 
            FALSE, c(-12, 12), seq(-12, 12, 4), 
            "Difference between\nfeeling thermometers\ntoward women and men", 
            c("Women", "Men"))
title("Heterogeneous treatment effect by gender, using all respondents (Japan)", 
      cex.main = 1.5, outer = TRUE)
dev.off()

# Figure A10
cairo_pdf("Figure_A10.pdf", width = 4.8, height = 1.4, pointsize = 7)
par(mar = c(1, 2.5, 4, 0.5), lwd = 0.5)
layout(matrix(1:4, 1, 4, byrow = TRUE))
effect.plot(subset(combined.data, T2 == 0), 
            "intv.gender.gap", "T1", "Japan", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Gender pay gap intervention\n", 
            c("Korea", "Japan"))
effect.plot(subset(combined.data, T2 == 0), 
            "intv.gender.policy", "T1", "Japan", 
            FALSE, c(-0.5, 1), seq(-0.5, 1, 0.25), 
            "Equal Employment Act\n", 
            c("Korea", "Japan"))
effect.plot(subset(combined.data, T2 == 0), 
            "feeling.women", "T1", "Japan", 
            FALSE, c(-10, 10), seq(-10, 10, 5), 
            "Feeling twd. women\n", 
            c("Korea", "Japan"))
effect.plot(subset(combined.data, T2 == 0), 
            "feeling.gender.diff", "T1", "Japan", 
            FALSE, c(-12, 12), seq(-12, 12, 4), 
            "Diff. b/w feeling\n", 
            c("Korea", "Japan"))
dev.off()

#### coefficient tables ####
## South Korea
# Table A6 (a)
H1.gap.intv.gap.Korea <- lm_robust(intv.income.gap ~ (T1 + T2) * income.gap.undest, 
                                   Korea.data)
coef.table(H1.gap.intv.gap.Korea, FALSE, FALSE)

H1.gap.intv.gap.ctrl.Korea <- lm_robust(intv.income.gap ~ (T1 + T2) * income.gap.undest + 
                                          factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                          income.incl.NA + ideology + HS + BS + MS, 
                                        Korea.data)
coef.table(H1.gap.intv.gap.ctrl.Korea, FALSE, FALSE)

H1.rank.intv.gap.Korea <- lm_robust(intv.income.gap ~ (T1 + T2) * income.rank.undest, 
                                    Korea.data)
coef.table(H1.rank.intv.gap.Korea, FALSE, FALSE)

H1.rank.intv.gap.ctrl.Korea <- lm_robust(intv.income.gap ~ (T1 + T2) * income.rank.undest + 
                                           factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                           income.incl.NA + ideology + HS + BS + MS, 
                                         Korea.data)
coef.table(H1.rank.intv.gap.ctrl.Korea, FALSE, FALSE)

# Table A6 (b)
H1.gap.intv.policy.Korea <- lm_robust(intv.income.policy ~ (T1 + T2) * income.gap.undest, 
                                      Korea.data)
coef.table(H1.gap.intv.policy.Korea, FALSE, FALSE)

H1.gap.intv.policy.ctrl.Korea <- lm_robust(intv.income.policy ~ (T1 + T2) * income.gap.undest + 
                                             factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS, 
                                           Korea.data)
coef.table(H1.gap.intv.policy.ctrl.Korea, FALSE, FALSE)

H1.rank.intv.policy.Korea <- lm_robust(intv.income.policy ~ (T1 + T2) * income.rank.undest, 
                                       Korea.data)
coef.table(H1.rank.intv.policy.Korea, FALSE, FALSE)

H1.rank.intv.policy.ctrl.Korea <- lm_robust(intv.income.policy ~ (T1 + T2) * income.rank.undest + 
                                              factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                              income.incl.NA + ideology + HS + BS + MS, 
                                            Korea.data)
coef.table(H1.rank.intv.policy.ctrl.Korea, FALSE, FALSE)

# Table A6 (c)
H1.gap.feeling.abs.Korea <- lm_robust(feeling.poor ~ (T1 + T2) * income.gap.undest, 
                                      Korea.data)
coef.table(H1.gap.feeling.abs.Korea, FALSE, FALSE)

H1.gap.feeling.abs.ctrl.Korea <- lm_robust(feeling.poor ~ (T1 + T2) * income.gap.undest + 
                                             factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS, 
                                           Korea.data)
coef.table(H1.gap.feeling.abs.ctrl.Korea, FALSE, FALSE)

H1.rank.feeling.abs.Korea <- lm_robust(feeling.poor ~ (T1 + T2) * income.rank.undest, 
                                       Korea.data)
coef.table(H1.rank.feeling.abs.Korea, FALSE, FALSE)

H1.rank.feeling.abs.ctrl.Korea <- lm_robust(feeling.poor ~ (T1 + T2) * income.rank.undest + 
                                              factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                              income.incl.NA + ideology + HS + BS + MS, 
                                            Korea.data)
coef.table(H1.rank.feeling.abs.ctrl.Korea, FALSE, FALSE)

# Table A6 (d)
H1.gap.feeling.rel.Korea <- lm_robust(feeling.income.diff ~ (T1 + T2) * income.gap.undest, 
                                      Korea.data)
coef.table(H1.gap.feeling.rel.Korea, FALSE, FALSE)

H1.gap.feeling.rel.ctrl.Korea <- lm_robust(feeling.income.diff ~ (T1 + T2) * income.gap.undest + 
                                             factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS, 
                                           Korea.data)
coef.table(H1.gap.feeling.rel.ctrl.Korea, FALSE, FALSE)

H1.rank.feeling.rel.Korea <- lm_robust(feeling.income.diff ~ (T1 + T2) * income.rank.undest, 
                                       Korea.data)
coef.table(H1.rank.feeling.rel.Korea, FALSE, FALSE)

H1.rank.feeling.rel.ctrl.Korea <- lm_robust(feeling.income.diff ~ (T1 + T2) * income.rank.undest + 
                                              factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                              income.incl.NA + ideology + HS + BS + MS, 
                                            Korea.data)
coef.table(H1.rank.feeling.rel.ctrl.Korea, FALSE, FALSE)

# Table A7 (a)
H2.gap.intv.gap.Korea <- lm_robust(intv.gender.gap ~ (T1 + T2) * gender.gap.undest, 
                                   Korea.data)
coef.table(H2.gap.intv.gap.Korea, FALSE, TRUE)

H2.gap.intv.gap.ctrl.Korea <- lm_robust(intv.gender.gap ~ (T1 + T2) * gender.gap.undest + 
                                          factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                          income.incl.NA + ideology + HS + BS + MS, 
                                        Korea.data)
coef.table(H2.gap.intv.gap.ctrl.Korea, FALSE, TRUE)

H2.rank.intv.gap.Korea <- lm_robust(intv.gender.gap ~ (T1 + T2) * gender.rank.undest, 
                                    Korea.data)
coef.table(H2.rank.intv.gap.Korea, FALSE, TRUE)

H2.rank.intv.gap.ctrl.Korea <- lm_robust(intv.gender.gap ~ (T1 + T2) * gender.rank.undest + 
                                           factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                           income.incl.NA + ideology + HS + BS + MS, 
                                         Korea.data)
coef.table(H2.rank.intv.gap.ctrl.Korea, FALSE, TRUE)

# Table A7 (b)
H2.gap.intv.policy.Korea <- lm_robust(intv.gender.policy ~ (T1 + T2) * gender.gap.undest, 
                                      Korea.data)
coef.table(H2.gap.intv.policy.Korea, FALSE, TRUE)

H2.gap.intv.policy.ctrl.Korea <- lm_robust(intv.gender.policy ~ (T1 + T2) * gender.gap.undest + 
                                             factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS, 
                                           Korea.data)
coef.table(H2.gap.intv.policy.ctrl.Korea, FALSE, TRUE)

H2.rank.intv.policy.Korea <- lm_robust(intv.gender.policy ~ (T1 + T2) * gender.rank.undest, 
                                       Korea.data)
coef.table(H2.rank.intv.policy.Korea, FALSE, TRUE)

H2.rank.intv.policy.ctrl.Korea <- lm_robust(intv.gender.policy ~ (T1 + T2) * gender.rank.undest + 
                                              factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                              income.incl.NA + ideology + HS + BS + MS, 
                                            Korea.data)
coef.table(H2.rank.intv.policy.ctrl.Korea, FALSE, TRUE)

# Table A7 (c)
H2.gap.feeling.abs.Korea <- lm_robust(feeling.women ~ (T1 + T2) * gender.gap.undest, 
                                      Korea.data)
coef.table(H2.gap.feeling.abs.Korea, FALSE, TRUE)

H2.gap.feeling.abs.ctrl.Korea <- lm_robust(feeling.women ~ (T1 + T2) * gender.gap.undest + 
                                             factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS, 
                                           Korea.data)
coef.table(H2.gap.feeling.abs.ctrl.Korea, FALSE, TRUE)

H2.rank.feeling.abs.Korea <- lm_robust(feeling.women ~ (T1 + T2) * gender.rank.undest, 
                                       Korea.data)
coef.table(H2.rank.feeling.abs.Korea, FALSE, TRUE)

H2.rank.feeling.abs.ctrl.Korea <- lm_robust(feeling.women ~ (T1 + T2) * gender.rank.undest + 
                                              factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                              income.incl.NA + ideology + HS + BS + MS, 
                                            Korea.data)
coef.table(H2.rank.feeling.abs.ctrl.Korea, FALSE, TRUE)

# Table A7 (d)
H2.gap.feeling.rel.Korea <- lm_robust(feeling.gender.diff ~ (T1 + T2) * gender.gap.undest, 
                                      Korea.data)
coef.table(H2.gap.feeling.rel.Korea, FALSE, TRUE)

H2.gap.feeling.rel.ctrl.Korea <- lm_robust(feeling.gender.diff ~ (T1 + T2) * gender.gap.undest + 
                                             factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS, 
                                           Korea.data)
coef.table(H2.gap.feeling.rel.ctrl.Korea, FALSE, TRUE)

H2.rank.feeling.rel.Korea <- lm_robust(feeling.gender.diff ~ (T1 + T2) * gender.rank.undest, 
                                       Korea.data)
coef.table(H2.rank.feeling.rel.Korea, FALSE, TRUE)

H2.rank.feeling.rel.ctrl.Korea <- lm_robust(feeling.gender.diff ~ (T1 + T2) * gender.rank.undest + 
                                              factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                              income.incl.NA + ideology + HS + BS + MS, 
                                            Korea.data)
coef.table(H2.rank.feeling.rel.ctrl.Korea, FALSE, TRUE)

# Table A8 (a)
H4.gap.intv.gap.Korea <- lm_robust(intv.income.gap ~ (T1 + T2) * high.income, 
                                   Korea.data, subset = income.gap.undest == 1)
coef.table(H4.gap.intv.gap.Korea, FALSE, FALSE)

H4.gap.intv.gap.ctrl.Korea <- lm_robust(intv.income.gap ~ (T1 + T2) * high.income + 
                                          factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                          income.incl.NA + ideology + HS + BS + MS, 
                                        Korea.data, subset = income.gap.undest == 1)
coef.table(H4.gap.intv.gap.ctrl.Korea, FALSE, FALSE)

H4.rank.intv.gap.Korea <- lm_robust(intv.income.gap ~ (T1 + T2) * high.income, 
                                    Korea.data, subset = income.rank.undest == 1)
coef.table(H4.rank.intv.gap.Korea, FALSE, FALSE)

H4.rank.intv.gap.ctrl.Korea <- lm_robust(intv.income.gap ~ (T1 + T2) * high.income + 
                                           factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                           income.incl.NA + ideology + HS + BS + MS, 
                                         Korea.data, subset = income.rank.undest == 1)
coef.table(H4.rank.intv.gap.ctrl.Korea, FALSE, FALSE)

# Table A8 (b)
H4.gap.intv.policy.Korea <- lm_robust(intv.income.policy ~ (T1 + T2) * high.income, 
                                      Korea.data, subset = income.gap.undest == 1)
coef.table(H4.gap.intv.policy.Korea, FALSE, FALSE)

H4.gap.intv.policy.ctrl.Korea <- lm_robust(intv.income.policy ~ (T1 + T2) * high.income + 
                                             factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS, 
                                           Korea.data, subset = income.gap.undest == 1)
coef.table(H4.gap.intv.policy.ctrl.Korea, FALSE, FALSE)

H4.rank.intv.policy.Korea <- lm_robust(intv.income.policy ~ (T1 + T2) * high.income, 
                                       Korea.data, subset = income.rank.undest == 1)
coef.table(H4.rank.intv.policy.Korea, FALSE, FALSE)

H4.rank.intv.policy.ctrl.Korea <- lm_robust(intv.income.policy ~ (T1 + T2) * high.income + 
                                              factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                              income.incl.NA + ideology + HS + BS + MS, 
                                            Korea.data, subset = income.rank.undest == 1)
coef.table(H4.rank.intv.policy.ctrl.Korea, FALSE, FALSE)

# Table A8 (c)
H4.gap.feeling.abs.Korea <- lm_robust(feeling.poor ~ (T1 + T2) * high.income, 
                                      Korea.data, subset = income.gap.undest == 1)
coef.table(H4.gap.feeling.abs.Korea, FALSE, FALSE)

H4.gap.feeling.abs.ctrl.Korea <- lm_robust(feeling.poor ~ (T1 + T2) * high.income + 
                                             factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS, 
                                           Korea.data, subset = income.gap.undest == 1)
coef.table(H4.gap.feeling.abs.ctrl.Korea, FALSE, FALSE)

H4.rank.feeling.abs.Korea <- lm_robust(feeling.poor ~ (T1 + T2) * high.income, 
                                       Korea.data, subset = income.rank.undest == 1)
coef.table(H4.rank.feeling.abs.Korea, FALSE, FALSE)

H4.rank.feeling.abs.ctrl.Korea <- lm_robust(feeling.poor ~ (T1 + T2) * high.income + 
                                              factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                              income.incl.NA + ideology + HS + BS + MS, 
                                            Korea.data, subset = income.rank.undest == 1)
coef.table(H4.rank.feeling.abs.ctrl.Korea, FALSE, FALSE)

# Table A8 (d)
H4.gap.feeling.rel.Korea <- lm_robust(feeling.income.diff ~ (T1 + T2) * high.income, 
                                      Korea.data, subset = income.gap.undest == 1)
coef.table(H4.gap.feeling.rel.Korea, FALSE, FALSE)

H4.gap.feeling.rel.ctrl.Korea <- lm_robust(feeling.income.diff ~ (T1 + T2) * high.income + 
                                             factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS, 
                                           Korea.data, subset = income.gap.undest == 1)
coef.table(H4.gap.feeling.rel.ctrl.Korea, FALSE, FALSE)

H4.rank.feeling.rel.Korea <- lm_robust(feeling.income.diff ~ (T1 + T2) * high.income, 
                                       Korea.data, subset = income.rank.undest == 1)
coef.table(H4.rank.feeling.rel.Korea, FALSE, FALSE)

H4.rank.feeling.rel.ctrl.Korea <- lm_robust(feeling.income.diff ~ (T1 + T2) * high.income + 
                                              factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                              income.incl.NA + ideology + HS + BS + MS, 
                                            Korea.data, subset = income.rank.undest == 1)
coef.table(H4.rank.feeling.rel.ctrl.Korea, FALSE, FALSE)

# Table A9 (a)
H5.gap.intv.gap.Korea <- lm_robust(intv.gender.gap ~ (T1 + T2) * male, 
                                   Korea.data, subset = gender.gap.undest == 1)
coef.table(H5.gap.intv.gap.Korea, FALSE, TRUE)

H5.gap.intv.gap.ctrl.Korea <- lm_robust(intv.gender.gap ~ (T1 + T2) * male + 
                                          age + I(age ^ 2) + factor(educ.cat) + 
                                          income.incl.NA + ideology + HS + BS + MS, 
                                        Korea.data, subset = gender.gap.undest == 1)
coef.table(H5.gap.intv.gap.ctrl.Korea, FALSE, TRUE)

H5.rank.intv.gap.Korea <- lm_robust(intv.gender.gap ~ (T1 + T2) * male, 
                                    Korea.data, subset = gender.rank.undest == 1)
coef.table(H5.rank.intv.gap.Korea, FALSE, TRUE)

H5.rank.intv.gap.ctrl.Korea <- lm_robust(intv.gender.gap ~ (T1 + T2) * male + 
                                           age + I(age ^ 2) + factor(educ.cat) + 
                                           income.incl.NA + ideology + HS + BS + MS, 
                                         Korea.data, subset = gender.rank.undest == 1)
coef.table(H5.rank.intv.gap.ctrl.Korea, FALSE, TRUE)

# Table A9 (b)
H5.gap.intv.policy.Korea <- lm_robust(intv.gender.policy ~ (T1 + T2) * male, 
                                      Korea.data, subset = gender.gap.undest == 1)
coef.table(H5.gap.intv.policy.Korea, FALSE, TRUE)

H5.gap.intv.policy.ctrl.Korea <- lm_robust(intv.gender.policy ~ (T1 + T2) * male + 
                                             age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS, 
                                           Korea.data, subset = gender.gap.undest == 1)
coef.table(H5.gap.intv.policy.ctrl.Korea, FALSE, TRUE)

H5.rank.intv.policy.Korea <- lm_robust(intv.gender.policy ~ (T1 + T2) * male, 
                                       Korea.data, subset = gender.rank.undest == 1)
coef.table(H5.rank.intv.policy.Korea, FALSE, TRUE)

H5.rank.intv.policy.ctrl.Korea <- lm_robust(intv.gender.policy ~ (T1 + T2) * male + 
                                              age + I(age ^ 2) + factor(educ.cat) + 
                                              income.incl.NA + ideology + HS + BS + MS, 
                                            Korea.data, subset = gender.rank.undest == 1)
coef.table(H5.rank.intv.policy.ctrl.Korea, FALSE, TRUE)

# Table A9 (c)
H5.gap.feeling.abs.Korea <- lm_robust(feeling.women ~ (T1 + T2) * male, 
                                      Korea.data, subset = gender.gap.undest == 1)
coef.table(H5.gap.feeling.abs.Korea, FALSE, TRUE)

H5.gap.feeling.abs.ctrl.Korea <- lm_robust(feeling.women ~ (T1 + T2) * male + 
                                             age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS, 
                                           Korea.data, subset = gender.gap.undest == 1)
coef.table(H5.gap.feeling.abs.ctrl.Korea, FALSE, TRUE)

H5.rank.feeling.abs.Korea <- lm_robust(feeling.women ~ (T1 + T2) * male, 
                                       Korea.data, subset = gender.rank.undest == 1)
coef.table(H5.rank.feeling.abs.Korea, FALSE, TRUE)

H5.rank.feeling.abs.ctrl.Korea <- lm_robust(feeling.women ~ (T1 + T2) * male + 
                                              age + I(age ^ 2) + factor(educ.cat) + 
                                              income.incl.NA + ideology + HS + BS + MS, 
                                            Korea.data, subset = gender.rank.undest == 1)
coef.table(H5.rank.feeling.abs.ctrl.Korea, FALSE, TRUE)

# Table A9 (d)
H5.gap.feeling.rel.Korea <- lm_robust(feeling.gender.diff ~ (T1 + T2) * male, 
                                      Korea.data, subset = gender.gap.undest == 1)
coef.table(H5.gap.feeling.rel.Korea, FALSE, TRUE)

H5.gap.feeling.rel.ctrl.Korea <- lm_robust(feeling.gender.diff ~ (T1 + T2) * male + 
                                             age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS, 
                                           Korea.data, subset = gender.gap.undest == 1)
coef.table(H5.gap.feeling.rel.ctrl.Korea, FALSE, TRUE)

H5.rank.feeling.rel.Korea <- lm_robust(feeling.gender.diff ~ (T1 + T2) * male, 
                                       Korea.data, subset = gender.rank.undest == 1)
coef.table(H5.rank.feeling.rel.Korea, FALSE, TRUE)

H5.rank.feeling.rel.ctrl.Korea <- lm_robust(feeling.gender.diff ~ (T1 + T2) * male + 
                                              age + I(age ^ 2) + factor(educ.cat) + 
                                              income.incl.NA + ideology + HS + BS + MS, 
                                            Korea.data, subset = gender.rank.undest == 1)
coef.table(H5.rank.feeling.rel.ctrl.Korea, FALSE, TRUE)

## Japan
# Table A10 (a)
H1.gap.intv.gap.Japan <- lm_robust(intv.income.gap ~ (T1 + T2) * income.gap.undest, 
                                   Japan.data, fixed_effects = block)
coef.table(H1.gap.intv.gap.Japan, TRUE, FALSE)

H1.gap.intv.gap.ctrl.Japan <- lm_robust(intv.income.gap ~ (T1 + T2) * income.gap.undest + 
                                          age + I(age ^ 2) + factor(educ.cat) + 
                                          income.incl.NA + ideology + HS + BS + MS, 
                                        Japan.data, fixed_effects = block)
coef.table(H1.gap.intv.gap.ctrl.Japan, TRUE, FALSE)

H1.rank.intv.gap.Japan <- lm_robust(intv.income.gap ~ (T1 + T2) * income.rank.undest, 
                                    Japan.data, fixed_effects = block)
coef.table(H1.rank.intv.gap.Japan, TRUE, FALSE)

H1.rank.intv.gap.ctrl.Japan <- lm_robust(intv.income.gap ~ (T1 + T2) * income.rank.undest + 
                                           age + I(age ^ 2) + factor(educ.cat) + 
                                           income.incl.NA + ideology + HS + BS + MS, 
                                         Japan.data, fixed_effects = block)
coef.table(H1.rank.intv.gap.ctrl.Japan, TRUE, FALSE)

# Table A10 (b)
H1.gap.intv.policy.Japan <- lm_robust(intv.income.policy ~ (T1 + T2) * income.gap.undest, 
                                      Japan.data, fixed_effects = block)
coef.table(H1.gap.intv.policy.Japan, TRUE, FALSE)

H1.gap.intv.policy.ctrl.Japan <- lm_robust(intv.income.policy ~ (T1 + T2) * income.gap.undest + 
                                             age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS, 
                                           Japan.data, fixed_effects = block)
coef.table(H1.gap.intv.policy.ctrl.Japan, TRUE, FALSE)

H1.rank.intv.policy.Japan <- lm_robust(intv.income.policy ~ (T1 + T2) * income.rank.undest, 
                                       Japan.data, fixed_effects = block)
coef.table(H1.rank.intv.policy.Japan, TRUE, FALSE)

H1.rank.intv.policy.ctrl.Japan <- lm_robust(intv.income.policy ~ (T1 + T2) * income.rank.undest + 
                                              age + I(age ^ 2) + factor(educ.cat) + 
                                              income.incl.NA + ideology + HS + BS + MS, 
                                            Japan.data, fixed_effects = block)
coef.table(H1.rank.intv.policy.ctrl.Japan, TRUE, FALSE)

# Table A10 (c)
H1.gap.feeling.abs.Japan <- lm_robust(feeling.poor ~ (T1 + T2) * income.gap.undest, 
                                      Japan.data, fixed_effects = block)
coef.table(H1.gap.feeling.abs.Japan, TRUE, FALSE)

H1.gap.feeling.abs.ctrl.Japan <- lm_robust(feeling.poor ~ (T1 + T2) * income.gap.undest + 
                                             age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS, 
                                           Japan.data, fixed_effects = block)
coef.table(H1.gap.feeling.abs.ctrl.Japan, TRUE, FALSE)

H1.rank.feeling.abs.Japan <- lm_robust(feeling.poor ~ (T1 + T2) * income.rank.undest, 
                                       Japan.data, fixed_effects = block)
coef.table(H1.rank.feeling.abs.Japan, TRUE, FALSE)

H1.rank.feeling.abs.ctrl.Japan <- lm_robust(feeling.poor ~ (T1 + T2) * income.rank.undest + 
                                              age + I(age ^ 2) + factor(educ.cat) + 
                                              income.incl.NA + ideology + HS + BS + MS, 
                                            Japan.data, fixed_effects = block)
coef.table(H1.rank.feeling.abs.ctrl.Japan, TRUE, FALSE)

# Table A10 (d)
H1.gap.feeling.rel.Japan <- lm_robust(feeling.income.diff ~ (T1 + T2) * income.gap.undest, 
                                      Japan.data, fixed_effects = block)
coef.table(H1.gap.feeling.rel.Japan, TRUE, FALSE)

H1.gap.feeling.rel.ctrl.Japan <- lm_robust(feeling.income.diff ~ (T1 + T2) * income.gap.undest + 
                                             age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS, 
                                           Japan.data, fixed_effects = block)
coef.table(H1.gap.feeling.rel.ctrl.Japan, TRUE, FALSE)

H1.rank.feeling.rel.Japan <- lm_robust(feeling.income.diff ~ (T1 + T2) * income.rank.undest, 
                                       Japan.data, fixed_effects = block)
coef.table(H1.rank.feeling.rel.Japan, TRUE, FALSE)

H1.rank.feeling.rel.ctrl.Japan <- lm_robust(feeling.income.diff ~ (T1 + T2) * income.rank.undest + 
                                              age + I(age ^ 2) + factor(educ.cat) + 
                                              income.incl.NA + ideology + HS + BS + MS, 
                                            Japan.data, fixed_effects = block)
coef.table(H1.rank.feeling.rel.ctrl.Japan, TRUE, FALSE)

# Table A11 (a)
H2.gap.intv.gap.Japan <- lm_robust(intv.gender.gap ~ (T1 + T2) * gender.gap.undest, 
                                   Japan.data, fixed_effects = block)
coef.table(H2.gap.intv.gap.Japan, TRUE, TRUE)

H2.gap.intv.gap.ctrl.Japan <- lm_robust(intv.gender.gap ~ (T1 + T2) * gender.gap.undest + 
                                          age + I(age ^ 2) + factor(educ.cat) + 
                                          income.incl.NA + ideology + HS + BS + MS, 
                                        Japan.data, fixed_effects = block)
coef.table(H2.gap.intv.gap.ctrl.Japan, TRUE, TRUE)

H2.rank.intv.gap.Japan <- lm_robust(intv.gender.gap ~ (T1 + T2) * gender.rank.undest, 
                                    Japan.data, fixed_effects = block)
coef.table(H2.rank.intv.gap.Japan, TRUE, TRUE)

H2.rank.intv.gap.ctrl.Japan <- lm_robust(intv.gender.gap ~ (T1 + T2) * gender.rank.undest + 
                                           age + I(age ^ 2) + factor(educ.cat) + 
                                           income.incl.NA + ideology + HS + BS + MS, 
                                         Japan.data, fixed_effects = block)
coef.table(H2.rank.intv.gap.ctrl.Japan, TRUE, TRUE)

# Table A11 (b)
H2.gap.intv.policy.Japan <- lm_robust(intv.gender.policy ~ (T1 + T2) * gender.gap.undest, 
                                      Japan.data, fixed_effects = block)
coef.table(H2.gap.intv.policy.Japan, TRUE, TRUE)

H2.gap.intv.policy.ctrl.Japan <- lm_robust(intv.gender.policy ~ (T1 + T2) * gender.gap.undest + 
                                             age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS, 
                                           Japan.data, fixed_effects = block)
coef.table(H2.gap.intv.policy.ctrl.Japan, TRUE, TRUE)

H2.rank.intv.policy.Japan <- lm_robust(intv.gender.policy ~ (T1 + T2) * gender.rank.undest, 
                                       Japan.data, fixed_effects = block)
coef.table(H2.rank.intv.policy.Japan, TRUE, TRUE)

H2.rank.intv.policy.ctrl.Japan <- lm_robust(intv.gender.policy ~ (T1 + T2) * gender.rank.undest + 
                                              age + I(age ^ 2) + factor(educ.cat) + 
                                              income.incl.NA + ideology + HS + BS + MS, 
                                            Japan.data, fixed_effects = block)
coef.table(H2.rank.intv.policy.ctrl.Japan, TRUE, TRUE)

# Table A11 (c)
H2.gap.feeling.abs.Japan <- lm_robust(feeling.women ~ (T1 + T2) * gender.gap.undest, 
                                      Japan.data, fixed_effects = block)
coef.table(H2.gap.feeling.abs.Japan, TRUE, TRUE)

H2.gap.feeling.abs.ctrl.Japan <- lm_robust(feeling.women ~ (T1 + T2) * gender.gap.undest + 
                                             age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS, 
                                           Japan.data, fixed_effects = block)
coef.table(H2.gap.feeling.abs.ctrl.Japan, TRUE, TRUE)

H2.rank.feeling.abs.Japan <- lm_robust(feeling.women ~ (T1 + T2) * gender.rank.undest, 
                                       Japan.data, fixed_effects = block)
coef.table(H2.rank.feeling.abs.Japan, TRUE, TRUE)

H2.rank.feeling.abs.ctrl.Japan <- lm_robust(feeling.women ~ (T1 + T2) * gender.rank.undest + 
                                              age + I(age ^ 2) + factor(educ.cat) + 
                                              income.incl.NA + ideology + HS + BS + MS, 
                                            Japan.data, fixed_effects = block)
coef.table(H2.rank.feeling.abs.ctrl.Japan, TRUE, TRUE)

# Table A11 (d)
H2.gap.feeling.rel.Japan <- lm_robust(feeling.gender.diff ~ (T1 + T2) * gender.gap.undest, 
                                      Japan.data, fixed_effects = block)
coef.table(H2.gap.feeling.rel.Japan, TRUE, TRUE)

H2.gap.feeling.rel.ctrl.Japan <- lm_robust(feeling.gender.diff ~ (T1 + T2) * gender.gap.undest + 
                                             age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS, 
                                           Japan.data, fixed_effects = block)
coef.table(H2.gap.feeling.rel.ctrl.Japan, TRUE, TRUE)

H2.rank.feeling.rel.Japan <- lm_robust(feeling.gender.diff ~ (T1 + T2) * gender.rank.undest, 
                                       Japan.data, fixed_effects = block)
coef.table(H2.rank.feeling.rel.Japan, TRUE, TRUE)

H2.rank.feeling.rel.ctrl.Japan <- lm_robust(feeling.gender.diff ~ (T1 + T2) * gender.rank.undest + 
                                              age + I(age ^ 2) + factor(educ.cat) + 
                                              income.incl.NA + ideology + HS + BS + MS, 
                                            Japan.data, fixed_effects = block)
coef.table(H2.rank.feeling.rel.ctrl.Japan, TRUE, TRUE)

# Table A12 (a)
H4.gap.intv.gap.Japan <- lm_robust(intv.income.gap ~ (T1 + T2) * high.income, 
                                   Japan.data, fixed_effects = block, 
                                   subset = income.gap.undest == 1)
coef.table(H4.gap.intv.gap.Japan, TRUE, FALSE)

H4.gap.intv.gap.ctrl.Japan <- lm_robust(intv.income.gap ~ (T1 + T2) * high.income + 
                                          age + I(age ^ 2) + factor(educ.cat) + 
                                          income.incl.NA + ideology + HS + BS + MS, 
                                        Japan.data, fixed_effects = block, 
                                        subset = income.gap.undest == 1)
coef.table(H4.gap.intv.gap.ctrl.Japan, TRUE, FALSE)

H4.rank.intv.gap.Japan <- lm_robust(intv.income.gap ~ (T1 + T2) * high.income, 
                                    Japan.data, fixed_effects = block, 
                                    subset = income.rank.undest == 1)
coef.table(H4.rank.intv.gap.Japan, TRUE, FALSE)

H4.rank.intv.gap.ctrl.Japan <- lm_robust(intv.income.gap ~ (T1 + T2) * high.income + 
                                           age + I(age ^ 2) + factor(educ.cat) + 
                                           income.incl.NA + ideology + HS + BS + MS, 
                                         Japan.data, fixed_effects = block, 
                                         subset = income.rank.undest == 1)
coef.table(H4.rank.intv.gap.ctrl.Japan, TRUE, FALSE)

# Table A12 (b)
H4.gap.intv.policy.Japan <- lm_robust(intv.income.policy ~ (T1 + T2) * high.income, 
                                      Japan.data, fixed_effects = block, 
                                      subset = income.gap.undest == 1)
coef.table(H4.gap.intv.policy.Japan, TRUE, FALSE)

H4.gap.intv.policy.ctrl.Japan <- lm_robust(intv.income.policy ~ (T1 + T2) * high.income + 
                                             age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS, 
                                           Japan.data, fixed_effects = block, 
                                           subset = income.gap.undest == 1)
coef.table(H4.gap.intv.policy.ctrl.Japan, TRUE, FALSE)

H4.rank.intv.policy.Japan <- lm_robust(intv.income.policy ~ (T1 + T2) * high.income, 
                                       Japan.data, fixed_effects = block, 
                                       subset = income.rank.undest == 1)
coef.table(H4.rank.intv.policy.Japan, TRUE, FALSE)

H4.rank.intv.policy.ctrl.Japan <- lm_robust(intv.income.policy ~ (T1 + T2) * high.income + 
                                              age + I(age ^ 2) + factor(educ.cat) + 
                                              income.incl.NA + ideology + HS + BS + MS, 
                                            Japan.data, fixed_effects = block, 
                                            subset = income.rank.undest == 1)
coef.table(H4.rank.intv.policy.ctrl.Japan, TRUE, FALSE)

# Table A12 (c)
H4.gap.feeling.abs.Japan <- lm_robust(feeling.poor ~ (T1 + T2) * high.income, 
                                      Japan.data, fixed_effects = block, 
                                      subset = income.gap.undest == 1)
coef.table(H4.gap.feeling.abs.Japan, TRUE, FALSE)

H4.gap.feeling.abs.ctrl.Japan <- lm_robust(feeling.poor ~ (T1 + T2) * high.income + 
                                             age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS, 
                                           Japan.data, fixed_effects = block, 
                                           subset = income.gap.undest == 1)
coef.table(H4.gap.feeling.abs.ctrl.Japan, TRUE, FALSE)

H4.rank.feeling.abs.Japan <- lm_robust(feeling.poor ~ (T1 + T2) * high.income, 
                                       Japan.data, fixed_effects = block, 
                                       subset = income.rank.undest == 1)
coef.table(H4.rank.feeling.abs.Japan, TRUE, FALSE)

H4.rank.feeling.abs.ctrl.Japan <- lm_robust(feeling.poor ~ (T1 + T2) * high.income + 
                                              age + I(age ^ 2) + factor(educ.cat) + 
                                              income.incl.NA + ideology + HS + BS + MS, 
                                            Japan.data, fixed_effects = block, 
                                            subset = income.rank.undest == 1)
coef.table(H4.rank.feeling.abs.ctrl.Japan, TRUE, FALSE)

# Table A12 (d)
H4.gap.feeling.rel.Japan <- lm_robust(feeling.income.diff ~ (T1 + T2) * high.income, 
                                      Japan.data, fixed_effects = block, 
                                      subset = income.gap.undest == 1)
coef.table(H4.gap.feeling.rel.Japan, TRUE, FALSE)

H4.gap.feeling.rel.ctrl.Japan <- lm_robust(feeling.income.diff ~ (T1 + T2) * high.income + 
                                             age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS, 
                                           Japan.data, fixed_effects = block, 
                                           subset = income.gap.undest == 1)
coef.table(H4.gap.feeling.rel.ctrl.Japan, TRUE, FALSE)

H4.rank.feeling.rel.Japan <- lm_robust(feeling.income.diff ~ (T1 + T2) * high.income, 
                                       Japan.data, fixed_effects = block, 
                                       subset = income.rank.undest == 1)
coef.table(H4.rank.feeling.rel.Japan, TRUE, FALSE)

H4.rank.feeling.rel.ctrl.Japan <- lm_robust(feeling.income.diff ~ (T1 + T2) * high.income + 
                                              age + I(age ^ 2) + factor(educ.cat) + 
                                              income.incl.NA + ideology + HS + BS + MS, 
                                            Japan.data, fixed_effects = block, 
                                            subset = income.rank.undest == 1)
coef.table(H4.rank.feeling.rel.ctrl.Japan, TRUE, FALSE)

# Table A13 (a)
H5.gap.intv.gap.Japan <- lm_robust(intv.gender.gap ~ (T1 + T2) * male, 
                                   Japan.data, fixed_effects = block, 
                                   subset = gender.gap.undest == 1)
coef.table(H5.gap.intv.gap.Japan, TRUE, TRUE)

H5.gap.intv.gap.ctrl.Japan <- lm_robust(intv.gender.gap ~ (T1 + T2) * male + 
                                          age + I(age ^ 2) + factor(educ.cat) + 
                                          income.incl.NA + ideology + HS + BS + MS, 
                                        Japan.data, fixed_effects = block, 
                                        subset = gender.gap.undest == 1)
coef.table(H5.gap.intv.gap.ctrl.Japan, TRUE, TRUE)

H5.rank.intv.gap.Japan <- lm_robust(intv.gender.gap ~ (T1 + T2) * male, 
                                    Japan.data, fixed_effects = block, 
                                    subset = gender.rank.undest == 1)
coef.table(H5.rank.intv.gap.Japan, TRUE, TRUE)

H5.rank.intv.gap.ctrl.Japan <- lm_robust(intv.gender.gap ~ (T1 + T2) * male + 
                                           age + I(age ^ 2) + factor(educ.cat) + 
                                           income.incl.NA + ideology + HS + BS + MS, 
                                         Japan.data, fixed_effects = block, 
                                         subset = gender.rank.undest == 1)
coef.table(H5.rank.intv.gap.ctrl.Japan, TRUE, TRUE)

# Table A13 (b)
H5.gap.intv.policy.Japan <- lm_robust(intv.gender.policy ~ (T1 + T2) * male, 
                                      Japan.data, fixed_effects = block, 
                                      subset = gender.gap.undest == 1)
coef.table(H5.gap.intv.policy.Japan, TRUE, TRUE)

H5.gap.intv.policy.ctrl.Japan <- lm_robust(intv.gender.policy ~ (T1 + T2) * male + 
                                             age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS, 
                                           Japan.data, fixed_effects = block, 
                                           subset = gender.gap.undest == 1)
coef.table(H5.gap.intv.policy.ctrl.Japan, TRUE, TRUE)

H5.rank.intv.policy.Japan <- lm_robust(intv.gender.policy ~ (T1 + T2) * male, 
                                       Japan.data, fixed_effects = block, 
                                       subset = gender.rank.undest == 1)
coef.table(H5.rank.intv.policy.Japan, TRUE, TRUE)

H5.rank.intv.policy.ctrl.Japan <- lm_robust(intv.gender.policy ~ (T1 + T2) * male + 
                                              age + I(age ^ 2) + factor(educ.cat) + 
                                              income.incl.NA + ideology + HS + BS + MS, 
                                            Japan.data, fixed_effects = block, 
                                            subset = gender.rank.undest == 1)
coef.table(H5.rank.intv.policy.ctrl.Japan, TRUE, TRUE)

# Table A13 (c)
H5.gap.feeling.abs.Japan <- lm_robust(feeling.women ~ (T1 + T2) * male, 
                                      Japan.data, fixed_effects = block, 
                                      subset = gender.gap.undest == 1)
coef.table(H5.gap.feeling.abs.Japan, TRUE, TRUE)

H5.gap.feeling.abs.ctrl.Japan <- lm_robust(feeling.women ~ (T1 + T2) * male + 
                                             age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS, 
                                           Japan.data, fixed_effects = block, 
                                           subset = gender.gap.undest == 1)
coef.table(H5.gap.feeling.abs.ctrl.Japan, TRUE, TRUE)

H5.rank.feeling.abs.Japan <- lm_robust(feeling.women ~ (T1 + T2) * male, 
                                       Japan.data, fixed_effects = block, 
                                       subset = gender.rank.undest == 1)
coef.table(H5.rank.feeling.abs.Japan, TRUE, TRUE)

H5.rank.feeling.abs.ctrl.Japan <- lm_robust(feeling.women ~ (T1 + T2) * male + 
                                              age + I(age ^ 2) + factor(educ.cat) + 
                                              income.incl.NA + ideology + HS + BS + MS, 
                                            Japan.data, fixed_effects = block, 
                                            subset = gender.rank.undest == 1)
coef.table(H5.rank.feeling.abs.ctrl.Japan, TRUE, TRUE)

# Table A13 (d)
H5.gap.feeling.rel.Japan <- lm_robust(feeling.gender.diff ~ (T1 + T2) * male, 
                                      Japan.data, fixed_effects = block, 
                                      subset = gender.gap.undest == 1)
coef.table(H5.gap.feeling.rel.Japan, TRUE, TRUE)

H5.gap.feeling.rel.ctrl.Japan <- lm_robust(feeling.gender.diff ~ (T1 + T2) * male + 
                                             age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS, 
                                           Japan.data, fixed_effects = block, 
                                           subset = gender.gap.undest == 1)
coef.table(H5.gap.feeling.rel.ctrl.Japan, TRUE, TRUE)

H5.rank.feeling.rel.Japan <- lm_robust(feeling.gender.diff ~ (T1 + T2) * male, 
                                       Japan.data, fixed_effects = block, 
                                       subset = gender.rank.undest == 1)
coef.table(H5.rank.feeling.rel.Japan, TRUE, TRUE)

H5.rank.feeling.rel.ctrl.Japan <- lm_robust(feeling.gender.diff ~ (T1 + T2) * male + 
                                              age + I(age ^ 2) + factor(educ.cat) + 
                                              income.incl.NA + ideology + HS + BS + MS, 
                                            Japan.data, fixed_effects = block, 
                                            subset = gender.rank.undest == 1)
coef.table(H5.rank.feeling.rel.ctrl.Japan, TRUE, TRUE)

## cross-country comparison
# Table A14 (a)
H6.gap.intv.gap <- lm_robust(intv.gender.gap ~ (T1 + T2) * Japan, 
                             combined.data, subset = gender.gap.undest == 1)
coef.table(H6.gap.intv.gap, FALSE, TRUE)

H6.gap.intv.gap.ctrl <- lm_robust(intv.gender.gap ~ (T1 + T2) * Japan + 
                                    factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                    income.incl.NA + ideology + HS + BS + MS, 
                                  combined.data, subset = gender.gap.undest == 1)
coef.table(H6.gap.intv.gap.ctrl, FALSE, TRUE)

H6.rank.intv.gap <- lm_robust(intv.gender.gap ~ (T1 + T2) * Japan, 
                              combined.data, subset = gender.rank.undest == 1)
coef.table(H6.rank.intv.gap, FALSE, TRUE)

H6.rank.intv.gap.ctrl <- lm_robust(intv.gender.gap ~ (T1 + T2) * Japan + 
                                     factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                     income.incl.NA + ideology + HS + BS + MS, 
                                   combined.data, subset = gender.rank.undest == 1)
coef.table(H6.rank.intv.gap.ctrl, FALSE, TRUE)

# Table A14 (b)
H6.gap.intv.policy <- lm_robust(intv.gender.policy ~ (T1 + T2) * Japan, 
                                combined.data, subset = gender.gap.undest == 1)
coef.table(H6.gap.intv.policy, FALSE, TRUE)

H6.gap.intv.policy.ctrl <- lm_robust(intv.gender.policy ~ (T1 + T2) * Japan + 
                                       factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                       income.incl.NA + ideology + HS + BS + MS, 
                                     combined.data, subset = gender.gap.undest == 1)
coef.table(H6.gap.intv.policy.ctrl, FALSE, TRUE)

H6.rank.intv.policy <- lm_robust(intv.gender.policy ~ (T1 + T2) * Japan, 
                                 combined.data, subset = gender.rank.undest == 1)
coef.table(H6.rank.intv.policy, FALSE, TRUE)

H6.rank.intv.policy.ctrl <- lm_robust(intv.gender.policy ~ (T1 + T2) * Japan + 
                                        factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                        income.incl.NA + ideology + HS + BS + MS, 
                                      combined.data, subset = gender.rank.undest == 1)
coef.table(H6.rank.intv.policy.ctrl, FALSE, TRUE)

# Table A14 (c)
H6.gap.feeling.abs <- lm_robust(feeling.women ~ (T1 + T2) * Japan, 
                                combined.data, subset = gender.gap.undest == 1)
coef.table(H6.gap.feeling.abs, FALSE, TRUE)

H6.gap.feeling.abs.ctrl <- lm_robust(feeling.women ~ (T1 + T2) * Japan + 
                                       factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                       income.incl.NA + ideology + HS + BS + MS, 
                                     combined.data, subset = gender.gap.undest == 1)
coef.table(H6.gap.feeling.abs.ctrl, FALSE, TRUE)

H6.rank.feeling.abs <- lm_robust(feeling.women ~ (T1 + T2) * Japan, 
                                 combined.data, subset = gender.rank.undest == 1)
coef.table(H6.rank.feeling.abs, FALSE, TRUE)

H6.rank.feeling.abs.ctrl <- lm_robust(feeling.women ~ (T1 + T2) * Japan + 
                                        factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                        income.incl.NA + ideology + HS + BS + MS, 
                                      combined.data, subset = gender.rank.undest == 1)
coef.table(H6.rank.feeling.abs.ctrl, FALSE, TRUE)

# Table A14 (d)
H6.gap.feeling.rel <- lm_robust(feeling.gender.diff ~ (T1 + T2) * Japan, 
                                combined.data, subset = gender.gap.undest == 1)
coef.table(H6.gap.feeling.rel, FALSE, TRUE)

H6.gap.feeling.rel.ctrl <- lm_robust(feeling.gender.diff ~ (T1 + T2) * Japan + 
                                       factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                       income.incl.NA + ideology + HS + BS + MS, 
                                     combined.data, subset = gender.gap.undest == 1)
coef.table(H6.gap.feeling.rel.ctrl, FALSE, TRUE)

H6.rank.feeling.rel <- lm_robust(feeling.gender.diff ~ (T1 + T2) * Japan, 
                                 combined.data, subset = gender.rank.undest == 1)
coef.table(H6.rank.feeling.rel, FALSE, TRUE)

H6.rank.feeling.rel.ctrl <- lm_robust(feeling.gender.diff ~ (T1 + T2) * Japan + 
                                        factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                        income.incl.NA + ideology + HS + BS + MS, 
                                      combined.data, subset = gender.rank.undest == 1)
coef.table(H6.rank.feeling.rel.ctrl, FALSE, TRUE)

#### cross-country comparison controlling for the interaction terms of covariates ####
# Table A15 (a)
H6.gap.intv.gap.ctrl.2 <- lm_robust(intv.gender.gap ~ T1 + T2 + Japan + 
                                      (T1 + T2) * 
                                      (factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                         income.incl.NA + ideology + HS + BS + MS) + 
                                      (T1 + T2) : Japan, 
                                    combined.data, subset = gender.gap.undest == 1)
coef.table(H6.gap.intv.gap.ctrl.2, FALSE, TRUE)

H6.rank.intv.gap.ctrl.2 <- lm_robust(intv.gender.gap ~ T1 + T2 + Japan + 
                                       (T1 + T2) * 
                                       (factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                          income.incl.NA + ideology + HS + BS + MS) + 
                                       (T1 + T2) : Japan, 
                                     combined.data, subset = gender.rank.undest == 1)
coef.table(H6.rank.intv.gap.ctrl.2, FALSE, TRUE)

H6.gap.intv.policy.ctrl.2 <- lm_robust(intv.gender.policy ~ T1 + T2 + Japan + 
                                         (T1 + T2) * 
                                         (factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                            income.incl.NA + ideology + HS + BS + MS) + 
                                         (T1 + T2) : Japan, 
                                       combined.data, subset = gender.gap.undest == 1)
coef.table(H6.gap.intv.policy.ctrl.2, FALSE, TRUE)

H6.rank.intv.policy.ctrl.2 <- lm_robust(intv.gender.policy ~ T1 + T2 + Japan + 
                                          (T1 + T2) * 
                                          (factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS) + 
                                          (T1 + T2) : Japan, 
                                        combined.data, subset = gender.rank.undest == 1)
coef.table(H6.rank.intv.policy.ctrl.2, FALSE, TRUE)

# Table A15 (b)
H6.gap.feeling.abs.ctrl.2 <- lm_robust(feeling.women ~ T1 + T2 + Japan + 
                                         (T1 + T2) * 
                                         (factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                            income.incl.NA + ideology + HS + BS + MS) + 
                                         (T1 + T2) : Japan, 
                                       combined.data, subset = gender.gap.undest == 1)
coef.table(H6.gap.feeling.abs.ctrl.2, FALSE, TRUE)

H6.rank.feeling.abs.ctrl.2 <- lm_robust(feeling.women ~ T1 + T2 + Japan + 
                                          (T1 + T2) * 
                                          (factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS) + 
                                          (T1 + T2) : Japan, 
                                        combined.data, subset = gender.rank.undest == 1)
coef.table(H6.rank.feeling.abs.ctrl.2, FALSE, TRUE)

H6.gap.feeling.rel.ctrl.2 <- lm_robust(feeling.gender.diff ~ T1 + T2 + Japan + 
                                         (T1 + T2) * 
                                         (factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                            income.incl.NA + ideology + HS + BS + MS) + 
                                         (T1 + T2) : Japan, 
                                       combined.data, subset = gender.gap.undest == 1)
coef.table(H6.gap.feeling.rel.ctrl.2, FALSE, TRUE)

H6.rank.feeling.rel.ctrl.2 <- lm_robust(feeling.gender.diff ~ T1 + T2 + Japan + 
                                          (T1 + T2) * 
                                          (factor(gender) + age + I(age ^ 2) + factor(educ.cat) + 
                                             income.incl.NA + ideology + HS + BS + MS) + 
                                          (T1 + T2) : Japan, 
                                        combined.data, subset = gender.rank.undest == 1)
coef.table(H6.rank.feeling.rel.ctrl.2, FALSE, TRUE)

#### three-way interaction ####
# Figure A11
cairo_pdf("Figure_A11.pdf", width = 4.8, height = 1.4, pointsize = 7)
par(mar = c(2, 2.5, 3, 0.5), lwd = 0.5)
layout(matrix(1:4, 1, 4, byrow = TRUE))
effect.plot.2(subset(combined.data, T2 == 0 & gender.gap.undest == 1), 
              "intv.gender.gap", "T1", "male", "Japan", 
              c(-0.5, 1), seq(-0.5, 1, 0.25), 
              "Gender pay gap intervention\n(Only underest. for ratio)", 
              c("Women\n", "Men\n"))
text(0.85, 0.75, "Korea")
text(1.15, -0.45, "Japan")
effect.plot.2(subset(combined.data, T2 == 0 & gender.gap.undest == 1), 
              "intv.gender.policy", "T1", "male", "Japan", 
              c(-0.5, 1), seq(-0.5, 1, 0.25), 
              "Equal Employment Act\n(Only underest. for ratio)", 
              c("Women\n", "Men\n"))
effect.plot.2(subset(combined.data, T2 == 0 & gender.rank.undest == 1), 
              "intv.gender.gap", "T1", "male", "Japan", 
              c(-0.5, 1), seq(-0.5, 1, 0.25), 
              "Gender pay gap intervention\n(Only underest. for rank)", 
              c("Women\n", "Men\n"))
effect.plot.2(subset(combined.data, T2 == 0 & gender.rank.undest == 1), 
              "intv.gender.policy", "T1", "male", "Japan", 
              c(-0.5, 1), seq(-0.5, 1, 0.25), 
              "Equal Employment Act\n(Only underest. for rank)", 
              c("Women\n", "Men\n"))
dev.off()